Straw – WGCNA

0.1 Preprocessing:

Basic dada2 pipeline for 16S:
truncLen=c(200, 180)
trimLeft=c(19,20)
maxEE=c(3,4)
truncQ=2
pool=pseudo
database - SILVA 138.1

Delete all bad taxa with parameters below:
annosighted on the phylum level +
chloroplasts/mithondria +
low reads/richness(based on the plot results)

tree from SEPP-QIIME2-plagin:

qiime tools import
–input-path ref_filt.fasta
–output-path rep_seq.qza
–type ‘FeatureData[Sequence]’ qiime fragment-insertion sepp
–i-representative-sequences rep_seq.qza   –i-reference-database sepp-refs-silva-128.qza
–o-tree insertion-tree.qza
–o-placements insertion-placements.qza
–p-threads 40
unzip -p insertion-tree.qza /data/ > tree.nwk

ITS preprocessing:

dada2 with cutadapt trimming
maxEE=c(5,8)
truncQ = 8
minLen = 50
pool=“pseudo”
UNITE version - 10.05.2021

ML tree from IQ-TREE programm:

mafft –auto FASTA > FASTA_ALIGHNED
iqtree -s FASTA_ALIGHNED -nt 30 -B 1000

0.2 Import

Import phyloseq object and libraries
Dataset - straw decomposition time series - from factor Day - 10 levels

library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(ANCOMBC)
library(heatmaply)
library(compositions)
library(igraph)
library(WGCNA)
require(DESeq2)
require(phyloseq)

ps.f <- readRDS("psf2")

0.3 functions

#phyloseq object to ampvis2 object
#https://gist.github.com/KasperSkytte/8d0ca4206a66be7ff6d76fc4ab8e66c6

phyloseq_to_ampvis2 <- function(physeq) {
  #check object for class
  if(!any(class(physeq) %in% "phyloseq"))
    stop("physeq object must be of class \"phyloseq\"", call. = FALSE)
  
  #ampvis2 requires taxonomy and abundance table, phyloseq checks for the latter
  if(is.null(physeq@tax_table))
    stop("No taxonomy found in the phyloseq object and is required for ampvis2", call. = FALSE)
  
  #OTUs must be in rows, not columns
  if(phyloseq::taxa_are_rows(physeq))
    abund <- as.data.frame(phyloseq::otu_table(physeq)@.Data)
  else
    abund <- as.data.frame(t(phyloseq::otu_table(physeq)@.Data))
  
  #tax_table is assumed to have OTUs in rows too
  tax <- phyloseq::tax_table(physeq)@.Data
  
  #merge by rownames (OTUs)
  otutable <- merge(
    abund,
    tax,
    by = 0,
    all.x = TRUE,
    all.y = FALSE,
    sort = FALSE
  )
  colnames(otutable)[1] <- "OTU"
  
  #extract sample_data (metadata)
  if(!is.null(physeq@sam_data)) {
    metadata <- data.frame(
      phyloseq::sample_data(physeq),
      row.names = phyloseq::sample_names(physeq), 
      stringsAsFactors = FALSE, 
      check.names = FALSE
    )
    
    #check if any columns match exactly with rownames
    #if none matched assume row names are sample identifiers
    samplesCol <- unlist(lapply(metadata, function(x) {
      identical(x, rownames(metadata))}))
    
    if(any(samplesCol)) {
      #error if a column matched and it's not the first
      if(!samplesCol[[1]])
        stop("Sample ID's must be in the first column in the sample metadata, please reorder", call. = FALSE)
    } else {
      #assume rownames are sample identifiers, merge at the end with name "SampleID"
      if(any(colnames(metadata) %in% "SampleID"))
        stop("A column in the sample metadata is already named \"SampleID\" but does not seem to contain sample ID's", call. = FALSE)
      metadata$SampleID <- rownames(metadata)
      
      #reorder columns so SampleID is the first
      metadata <- metadata[, c(which(colnames(metadata) %in% "SampleID"), 1:(ncol(metadata)-1L)), drop = FALSE]
    }
  } else
    metadata <- NULL
  
  #extract phylogenetic tree, assumed to be of class "phylo"
  if(!is.null(physeq@phy_tree)) {
    tree <- phyloseq::phy_tree(physeq)
  } else
    tree <- NULL
  
  #extract OTU DNA sequences, assumed to be of class "XStringSet"
  if(!is.null(physeq@refseq)) {
    #convert XStringSet to DNAbin using a temporary file (easiest)
    fastaTempFile <- tempfile(pattern = "ampvis2_", fileext = ".fa")
    Biostrings::writeXStringSet(physeq@refseq, filepath = fastaTempFile)
  } else
    fastaTempFile <- NULL
  
  #load as normally with amp_load
  ampvis2::amp_load(
    otutable = otutable,
    metadata = metadata,
    tree = tree,
    fasta = fastaTempFile
  )
}

#variance stabilisation from DESeq2

ps_vst <- function(ps, factor){
  

  
  diagdds = phyloseq_to_deseq2(ps, as.formula(paste( "~", factor)))              
  diagdds = estimateSizeFactors(diagdds, type="poscounts")
  diagdds = estimateDispersions(diagdds, fitType = "local") 
  pst <- varianceStabilizingTransformation(diagdds)
  pst.dimmed <- t(as.matrix(assay(pst))) 
  # pst.dimmed[pst.dimmed < 0.0] <- 0.0
  ps.varstab <- ps
  otu_table(ps.varstab) <- otu_table(pst.dimmed, taxa_are_rows = FALSE) 
  return(ps.varstab)
}

#WGCNA visualisation
#result - list class object with attributes:
# ps - phyloseq object
# amp - ampwis2 object
# heat - heatmap with absalute read numbers
# heat_rel - hetmap with relative abundances
# tree - phylogenetic tree with taxonomy

color_filt <- function(ps, df){
  
      library(tidyverse)
      library(reshape2)
      library(gridExtra)

  
  l = list()
  for (i in levels(df$module)){
    message(i)
    color_name <-  df %>% 
      filter(module == i) %>% 
      pull(asv) %>% 
      unique()
    ps.col <- prune_taxa(color_name, ps)
    amp.col <- phyloseq_to_ampvis2(ps.col)
    heat <- amp_heatmap(amp.col, tax_show = 60, 
              group_by = "Day", 
              tax_aggregate = "OTU",
              tax_add = "Genus", 
              normalise=FALSE, 
              showRemainingTaxa = TRUE)
    
    ps.rel  <-  phyloseq::transform_sample_counts(ps, function(x) x / sum(x) * 100)
    ps.rel.col <- prune_taxa(color_name, ps.rel)
    amp.r <- phyloseq_to_ampvis2(ps.rel.col)
    heat.rel <- amp_heatmap(amp.r, tax_show = 60, 
          group_by = "Day", 
          tax_aggregate = "OTU",
          tax_add = "Genus", 
          normalise=FALSE, 
          showRemainingTaxa = TRUE)
    
    tree <- ps.col@phy_tree
    taxa <- as.data.frame(ps.col@tax_table@.Data)
    p1 <- ggtree(tree) + 
              geom_tiplab(size=2, align=TRUE, linesize=.5) +
              theme_tree2()
    
    taxa[taxa == "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"] <- "Allorhizobium"
    taxa[taxa == "Burkholderia-Caballeronia-Paraburkholderia"] <- "Burkholderia"
  
    tx <- taxa %>% 
      rownames_to_column("id") %>% 
      mutate(id = factor(id, levels = rev(get_taxa_name(p1)))) %>% 
      dplyr::select(-c(Kingdom, Species, Order)) %>% 
      melt(id.var = 'id')

    p2 <- ggplot(tx, aes(variable, id)) + 
      geom_tile(aes(fill = value), alpha = 0.4) +
      geom_text(aes(label = value), size = 3) +
      theme_bw() + 
      theme(legend.position = "none",
             axis.ticks.x=element_blank(),
             axis.text.y=element_blank(),
             axis.ticks.y=element_blank(),
             axis.title.x = element_blank(),
            axis.title.y = element_blank()) 

    p <- ggpubr::ggarrange(p1, p2, widths = c(0.9, 1))
    l[[i]] <- list("ps" = ps.col, 
                   "amp" = amp.col,
                   "heat" = heat,
                   "heat_rel" = heat.rel,
                   "tree" = p,
                   "taxa" = knitr::kable(taxa))
                    
  }
  return(l)
}


color_filt_broken <- function(ps, df, ps.pruned){
  
      library(tidyverse)
      library(reshape2)
      library(gridExtra)

  
  l = list()
  for (i in levels(df$module)){
    message(i)
    color_name <-  df %>% 
      filter(module == i) %>% 
      pull(asv) %>% 
      unique()
    ps.col <- prune_taxa(color_name, ps)
    ps.col.pruned <- prune_taxa(color_name, ps.pruned)
    amp.col <- phyloseq_to_ampvis2(ps.col)
    heat <- amp_heatmap(amp.col, tax_show = 60, 
              group_by = "Day", 
              tax_aggregate = "OTU",
              tax_add = "Genus", 
              normalise=FALSE, 
              showRemainingTaxa = TRUE)
    
    ps.rel  <-  phyloseq::transform_sample_counts(ps.col, function(x) x / sum(x) * 100)
    amp.r <- phyloseq_to_ampvis2(ps.rel)
    heat.rel <- amp_heatmap(amp.r, tax_show = 60, 
          group_by = "Day", 
          tax_aggregate = "OTU",
          tax_add = "Genus", 
          normalise=FALSE, 
          showRemainingTaxa = TRUE)
    
    tree <- ps.col.pruned@phy_tree
    taxa <- as.data.frame(ps.col.pruned@tax_table@.Data)
    p1 <- ggtree(tree) + 
              geom_tiplab(size=2, align=TRUE, linesize=.5) +
              theme_tree2()
    
    taxa[taxa == "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"] <- "Allorhizobium"
    taxa[taxa == "Burkholderia-Caballeronia-Paraburkholderia"] <- "Burkholderia"
  
    tx <- taxa %>% 
      rownames_to_column("id") %>% 
      mutate(id = factor(id, levels = rev(get_taxa_name(p1)))) %>% 
      dplyr::select(-c(Kingdom, Species, Order)) %>% 
      melt(id.var = 'id')

    p2 <- ggplot(tx, aes(variable, id)) + 
      geom_tile(aes(fill = value), alpha = 0.4) +
      geom_text(aes(label = value), size = 3) +
      theme_bw() + 
      theme(legend.position = "none",
             axis.ticks.x=element_blank(),
             axis.text.y=element_blank(),
             axis.ticks.y=element_blank(),
             axis.title.x = element_blank(),
            axis.title.y = element_blank()) 

    p <- ggpubr::ggarrange(p1, p2, widths = c(0.9, 1))
    l[[i]] <- list("ps" = ps.col, 
                   "amp" = amp.col,
                   "heat" = heat,
                   "heat_rel" = heat.rel,
                   "tree" = p,
                   "taxa" = knitr::kable(taxa))
                    
  }
  return(l)
}

detachAllPackages <- function() {

  basic.packages <- c("package:stats","package:graphics","package:grDevices","package:utils","package:datasets","package:methods","package:base")

  package.list <- search()[ifelse(unlist(gregexpr("package:",search()))==1,TRUE,FALSE)]

  package.list <- setdiff(package.list,basic.packages)

  if (length(package.list)>0)  for (package in package.list) detach(package, character.only=TRUE, force = TRUE)

}

plot_alpha_w_toc <- function(ps, group, metric) {
  
  require(phyloseq)
  require(ggplot2)
  
  ps_a <- prune_taxa(taxa_sums(ps) > 0, ps)
  
  er <- estimate_richness(ps_a)
  df_er <- cbind(ps_a@sam_data, er)
  df_er <- df_er %>% select(c(group, metric))
  stat.test <- aov(as.formula(paste0(metric, "~", group)), data = df_er) %>%
    rstatix::tukey_hsd()
  y <-  seq(max(er[[metric]]), length=length(stat.test$p.adj), by=max(er[[metric]]/20))

  plot_richness(ps_a, x=group, measures=metric) + 
    geom_boxplot() +
    geom_point(size=1.2, alpha=0.3) +
    ggpubr::stat_pvalue_manual(
      stat.test, 
      label = "p.adj.signif", 
      y.position = y) +
    theme_light() + 
    scale_color_brewer(palette="Dark2") +
    theme(axis.text.x = element_text(angle = 90),
          axis.title.x=element_blank()) +
    labs(y=paste(metric, "index")) 
}


# standart NMDS plot tool frop phyloseq with some additional aestatics 
# have stress value on plot - may work as fuck
beta_custom_norm_NMDS_elli_w <- function(ps, seed = 7888, normtype="vst", Color="What", Group="Repeat"){
  require(phyloseq)
  require(ggplot2)
  require(ggpubr)
  library(ggforce)
  
  
  ordination.b <- ordinate(ps, "NMDS", "bray")
  mds <- as.data.frame(ordination.b$points)
  p  <-  plot_ordination(ps,
                         ordination.b,
                         type="sample",
                         color = Color,
                         title="NMDS - Bray-Curtis",
                         # title=NULL,
                         axes = c(1,2) ) + 
    theme_bw() + 
    theme(text = element_text(size = 10)) + 
    geom_point(size = 3) +
    annotate("text",
    x=min(mds$MDS1) + abs(min(mds$MDS1))/7,
    y=max(mds$MDS2),
    label=paste0("Stress -- ", round(ordination.b$stress, 3))) +
    geom_mark_ellipse(aes_string(group = Group, label = Group),
                      label.fontsize = 10,
                      label.buffer = unit(2, "mm"),
                      label.minwidth = unit(5, "mm"),
                      con.cap = unit(0.1, "mm"),
                      con.colour='gray') +
    theme(legend.position = "none") +
    ggplot2::scale_fill_viridis_c(option = "H")
  
  return(p)
}


# alpha with aov + tukie post-hock - useless, but it looks pretty good
plot_alpha_w_toc <- function(ps, group, metric) {
  
  require(phyloseq)
  require(ggplot2)
  
  ps_a <- prune_taxa(taxa_sums(ps) > 0, ps)
  
  er <- estimate_richness(ps_a)
  df_er <- cbind(ps_a@sam_data, er)
  df_er <- df_er %>% select(c(group, metric))
  stat.test <- aov(as.formula(paste0(metric, "~", group)), data = df_er) %>%
    rstatix::tukey_hsd()
  y <- seq(max(er[[metric]]), length=length(stat.test$p.adj.signif[stat.test$p.adj.signif != "ns"]), by=max(er[[metric]]/20))
  
  plot_richness(ps_a, x=group, measures=metric) + 
    geom_boxplot() +
    geom_point(size=1.2, alpha=0.3) +
    stat_pvalue_manual(
      stat.test, 
      label = "p.adj.signif", 
      y.position = y,
      hide.ns=TRUE) +
    theme_light() + 
    scale_color_brewer(palette="Dark2") +
    theme(axis.text.x = element_text(angle = 90),
          axis.title.x=element_blank()) +
    labs(y=paste(metric, "index")) 
}

0.3.1 modify metadata

Add Group parameter to metadata -
- early - D01, D03, D05
- middle - D07, D08, D10
- late - D13, D14, D15

sample.data <- ps.f@sam_data %>% 
  data.frame() %>% 
  mutate(Group = if_else(Day %in% c("D01", "D03", "D05"), "early", 
                         if_else(Day %in% c("D07", "D08","D10"), "middle", "late"))) %>% 
  rename('Day' = 'Bag') %>% 
  mutate(Group = factor(Group, levels=c("early", "middle","late"))) %>% 
  mutate(Day = Bag %>% 
           forcats::fct_recode( "3" = "D01",
                                "14" = "D03",
                                "28" = "D05",
                                "49" = "D07" ,
                                "63" = "D08",
                                "91" = "D10",
                                "119" = "D12",
                                "140" = "D13",
                                "161" = "D14",
                                "182" = "D15")
         ) %>% 
  phyloseq::sample_data()

sample_data(ps.f) <- sample.data

0.4 Alpha diversity

p1 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "Observed")
p2 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "Shannon")
p3 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "InvSimpson")
ggarrange(p1, p2, p3, nrow = 1)

p.observed <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("Observed")) + 
  theme(axis.title.y = element_blank())

p.shannon <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("Shannon")) + 
  theme(axis.title.y = element_blank())

p.simpson <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("InvSimpson")) + 
  theme(axis.title.y = element_blank())

ggpubr::ggarrange(p.observed, p.shannon, p.simpson, ncol = 3)

0.4.0.1 mpd

mpd - index of alpha diversity “tree shagginess” - counted separately because it is in a separate package with built-in validity (based on permutations)
here is a link about the meanin of values in the columns:
https://www.rdocumentation.org/packages/picante/versions/1.8.2/topics/ses.mpd
in general the early stage is significantly less diverse than the other stages

physeq_merged <- merge_samples(ps.f, "Group", fun=sum)
# ps.f@sam_data

# picante::mpd(l_vst$blue$ps@otu_table@.Data, cophenetic(l_vst$blue$ps@phy_tree)) %>% 
#   mean(na.rm = TRUE)
# 
# picante::mpd(l_vst$salmon$ps@otu_table@.Data, cophenetic(l_vst$salmon$ps@phy_tree)) %>% 
#   mean(na.rm = TRUE)

mpd.res <- picante::ses.mpd(physeq_merged@otu_table@.Data, cophenetic(physeq_merged@phy_tree)) 

as.data.frame(mpd.res)

##        ntaxa  mpd.obs mpd.rand.mean mpd.rand.sd mpd.obs.rank mpd.obs.z
## early    364 1.742096      1.865821  0.02991525            1 -4.135846
## middle   490 1.839665      1.864686  0.02345011          151 -1.066981
## late     822 1.877129      1.864630  0.01354421          811  0.922841
##        mpd.obs.p runs
## early      0.001  999
## middle     0.151  999
## late       0.811  999

0.5 beta statistic

0.5.1 permanova

permanova - group significantlly different - dispersion between is more, than inside groups

dist <- phyloseq::distance(ps.f, "bray")
metadata <- as(sample_data(ps.f@sam_data), "data.frame")
vegan::adonis2(dist ~ Group, data = metadata)

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = dist ~ Group, data = metadata)
##          Df SumOfSqs      R2      F Pr(>F)    
## Group     2   3.2742 0.33893 8.2033  0.001 ***
## Residual 32   6.3861 0.66107                  
## Total    34   9.6604 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

0.5.2 dynamic of beta diversity

Bray Curtis distance by steps between each day.   (For example, all D15 versus D14.) Can we also put soil respiration on this picture as well?

ps.f.r <- rarefy_even_depth(ps.f, rngseed = 777)

avg.r <- ps.f.r@otu_table %>% 
  as.data.frame() %>%
  vegan::avgdist(10)

avg <- ps.f@otu_table %>% 
  as.data.frame() %>%
  vegan::avgdist(10)

# avg %>%
#   as.matrix() %>%
#   as_tibble(rownames= "sample") %>%
#   pivot_longer(-sample) %>%
#   filter(sample < name) %>%
#   mutate(repeat_a = str_replace(sample, ".*-", ""),
#          repeat_b = str_replace(name, ".*-", ""), 
#          day_a = as.numeric(str_replace(sapply(strsplit(sample, "-"), `[`, 3), "D", "")),
#           day_b = as.numeric(str_replace(sapply(strsplit(name, "-"), `[`, 3), "D", "")),
#          diff = abs(day_a - day_b),
#          early = day_a < 10) %>% 
#   filter(repeat_a == repeat_b & diff < 10) %>%
#   group_by(diff, repeat_a, early) %>%
#   summarize(median = median(value)) %>%
#   ungroup() %>%
#   ggplot(aes(x=diff, y=median, color=early, group=paste0(repeat_a, early))) +
#   geom_line(size=0.25) +
#   geom_smooth(aes(group=early), se=FALSE, size=4) +
#   labs(x="Distance between time points",
#        y="Median Bray-Curtis distance") +
#   scale_x_continuous(breaks=1:9) +
#   scale_color_manual(name=NULL,
#                      breaks=c(TRUE, FALSE),
#                      values=c("blue", "red"),
#                      labels=c("Early", "Late")) +
#   guides(color = guide_legend(override.aes = list(size=1))) +
#   theme_classic()

avg.r %>%
  as.matrix() %>%
  as_tibble(rownames= "sample") %>%
  pivot_longer(-sample) %>%
  filter(sample < name) %>%
  mutate(repeat_a = str_replace(sample, ".*-", ""),
         repeat_b = str_replace(name, ".*-", ""), 
         day_a = as.numeric(str_replace(sapply(strsplit(sample, "-"), `[`, 3), "D", "")),
          day_b = as.numeric(str_replace(sapply(strsplit(name, "-"), `[`, 3), "D", ""))) %>% 
  mutate(day_a = as.numeric(as.factor(day_a) %>% forcats::fct_recode("2" = "3", "3" = "5", "4" = "7", "5" = "8", "6" = "10", "7" = "13", "8" = "14", "9" = "14", "10" = "15")),
         day_b = as.numeric(as.factor(day_b) %>% forcats::fct_recode("2" = "3", "3" = "5", "4" = "7", "5" = "8", "6" = "10", "7" = "12", "8" = "13", "9" = "14", "10" = "15")),
         diff = abs(day_a - day_b)) %>% 
  filter(diff == 1) %>% 
  mutate(day_b = as.factor(day_b) %>% forcats::fct_recode("D03" = "2", "D05" = "3","D07" = "4",  "D08" = "5", "D10" = "6", "D13" = "7", "D14" = "8", "D15" = "10", "D12" = "7")) %>% 
  ggplot(aes(x=day_b, y=value)) +
    geom_boxplot() +
  theme_bw() +
  labs(x="Time points",
       y="Bray-Curtis distance")

0.5.3 basic NMDS plot

From the picture it looks like there is a gradual slowing down of changes -
but from the previous picture it follows that this is not the case - there is a dip between 7-8-10,
but on average the points differ quite consistently.

beta_custom_norm_NMDS_elli_w(ps.f, C="Group", G="Day")

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1048129 
## Run 1 stress 0.1372731 
## Run 2 stress 0.1318495 
## Run 3 stress 0.1048129 
## ... Procrustes: rmse 0.00004856956  max resid 0.0001949092 
## ... Similar to previous best
## Run 4 stress 0.1048129 
## ... New best solution
## ... Procrustes: rmse 0.00001842645  max resid 0.00006948497 
## ... Similar to previous best
## Run 5 stress 0.1426236 
## Run 6 stress 0.1048204 
## ... Procrustes: rmse 0.001770292  max resid 0.007527799 
## ... Similar to previous best
## Run 7 stress 0.1064532 
## Run 8 stress 0.1064733 
## Run 9 stress 0.1048204 
## ... Procrustes: rmse 0.001767644  max resid 0.007510715 
## ... Similar to previous best
## Run 10 stress 0.1064532 
## Run 11 stress 0.1318495 
## Run 12 stress 0.1064733 
## Run 13 stress 0.1064733 
## Run 14 stress 0.1460032 
## Run 15 stress 0.1048204 
## ... Procrustes: rmse 0.001769848  max resid 0.007525048 
## ... Similar to previous best
## Run 16 stress 0.1064532 
## Run 17 stress 0.1372728 
## Run 18 stress 0.1048129 
## ... Procrustes: rmse 0.000002711789  max resid 0.0000100797 
## ... Similar to previous best
## Run 19 stress 0.1460029 
## Run 20 stress 0.146003 
## *** Solution reached

0.6 Splitting the dataset

0.6.1 Basic information about dataset

amp <- phyloseq_to_ampvis2(ps.f)
amp

## ampvis2 object with 5 elements. 
## Summary of OTU table:
##      Samples         OTUs  Total#Reads    Min#Reads    Max#Reads Median#Reads 
##           35         1245       521643         4412        33561        12894 
##    Avg#Reads 
##     14904.09 
## 
## Assigned taxonomy:
##      Kingdom       Phylum        Class        Order       Family        Genus 
##   1245(100%)   1245(100%) 1233(99.04%)  1184(95.1%) 1063(85.38%)  736(59.12%) 
##      Species 
##    90(7.23%) 
## 
## Metadata variables: 5 
##  SampleID, Bag, Description, Group, Day

amp_heatmap(amp,
            tax_show = 22,
            group_by = "Day",
            tax_aggregate = "Phylum",
            tax_class = "Proteobacteriota",
            normalise=TRUE,
            showRemainingTaxa = TRUE)

## Warning: Transformation introduced infinite values in discrete y-axis

Let’s divide the dataset into two groups.

1st - more than 10% of samples should contain at least 10 reads - this group will be analyzed further

ps.inall <- phyloseq::filter_taxa(ps.f, function(x) sum(x > 10) > (0.1*length(x)), TRUE)

amp.inall <- phyloseq_to_ampvis2(ps.inall)
amp_heatmap(amp.inall,
            tax_show = 40,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=FALSE,
            showRemainingTaxa = TRUE)

The second part - the remaining, but more than 100 reeds in all samples(hereinafter - outliers majors(OM))
The remaining phylotypes are dropped from the analysis.

ps.exl  <- phyloseq::filter_taxa(ps.f, function(x) sum(x > 10) < (0.1*length(x)), TRUE)
ps.exl  <- prune_taxa(taxa_sums(ps.exl) > 100, ps.exl)
amp.exl <- phyloseq_to_ampvis2(ps.exl)
amp_heatmap(amp.exl,
            tax_show = 40,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=FALSE,
            showRemainingTaxa = TRUE)

Same, but with relative abundance.

ps.per  <-  phyloseq::transform_sample_counts(ps.f, function(x) x / sum(x) * 100) 
ps.exl.taxa <- taxa_names(ps.exl)
ps.per.exl <- prune_taxa(ps.exl.taxa, ps.per)
amp.exl.r <- phyloseq_to_ampvis2(ps.per.exl)

amp_heatmap(amp.exl.r, tax_show = 60, 
      group_by = "Day", 
      tax_aggregate = "OTU",
      tax_add = "Genus", 
      round =  2, 
      normalise=FALSE, 
      showRemainingTaxa = TRUE)

amp_heatmap(amp.exl.r, tax_show = 60, 
      group_by = "Day", 
      tax_aggregate = "OTU",
      tax_add = "Genus", 
      round =  2, 
      normalise=FALSE, )

OM aggregated by genus. Relative abundance.

amp_heatmap(amp.exl.r,
            tax_show = 30, 
            group_by = "Day", 
            tax_aggregate = "Genus",
            tax_add = "Phylum",
            tax_class = "Proteobacteria",
            round =  2,
            normalise=FALSE, 
            showRemainingTaxa = TRUE)

If we look at how OMs are distributed by day - it seems that the distribution of OMs is related not only to the specificity of individual bags, but also to purely technical features - the outlier values are mostly not with biological repeats, but with technical ones (red line - selected visually)

p_box <- phyloseq::sample_sums(ps.per.exl) %>% 
  as.data.frame(col.names = "values") %>% 
  setNames(., nm = "values") %>% 
  rownames_to_column("samples") %>% 
  mutate(Day = sapply(strsplit(samples, "-"), `[`, 3)) %>% 
  ggplot(aes(x=Day, y=values, color=Day, fill = Day)) +
  geom_boxplot(aes(color=Day, fill = Day)) +
  geom_point(color = "black", position = position_dodge(width=0.2)) +
  geom_hline(yintercept = 10, colour = "red") +
  theme_bw() +
  theme(legend.position = "none")
     
p_box <- p_box + viridis::scale_color_viridis(option = "H", discrete = TRUE, direction=1, begin=0.1, end = 0.9, alpha = 0.5)

p_box + viridis::scale_fill_viridis(option = "H", discrete = TRUE, direction=1, begin=0.1, end = 0.9, alpha = 0.3)

0.7 WGCNA

0.7.1 clr

after clr normalization / looks disgusting - let’s try vst normalization from DESeq2

otu.inall <- as.data.frame(ps.inall@otu_table@.Data) 
ps.inall.clr <- ps.inall
otu_table(ps.inall.clr) <- phyloseq::otu_table(compositions::clr(otu.inall), taxa_are_rows = FALSE)
data <- ps.inall.clr@otu_table@.Data %>% 
  as.data.frame()
rownames(data) <- as.character(ps.inall.clr@sam_data$Description)
powers <-  c(c(1:10), seq(from = 12, to=30, by=1))
sft <-  pickSoftThreshold(data, powerVector = powers, verbose = 5, networkType = "signed")

## pickSoftThreshold: will use block size 338.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 338 of 338

## Warning: executing %dopar% sequentially: no parallel backend registered

##    Power SFT.R.sq  slope truncated.R.sq  mean.k. median.k.  max.k.
## 1      1    0.185  23.60         0.8790 168.0000 168.00000 172.000
## 2      2    0.199 -15.50         0.8970  87.6000  87.20000  95.100
## 3      3    0.518 -13.10         0.8580  47.3000  46.80000  55.900
## 4      4    0.390  -6.46         0.8050  26.5000  26.00000  34.500
## 5      5    0.262  -3.46         0.7820  15.4000  14.90000  22.300
## 6      6    0.232  -2.30         0.7770   9.2300   8.78000  15.000
## 7      7    0.251  -1.87         0.6990   5.7300   5.36000  10.500
## 8      8    0.319  -1.75         0.7570   3.6800   3.41000   7.620
## 9      9    0.527  -2.07         0.7770   2.4400   2.20000   5.890
## 10    10    0.648  -2.28         0.8490   1.6600   1.47000   4.690
## 11    12    0.827  -2.41         0.8970   0.8420   0.69300   3.190
## 12    13    0.793  -2.38         0.7930   0.6220   0.49000   2.700
## 13    14    0.277  -3.97         0.0914   0.4690   0.35500   2.320
## 14    15    0.282  -3.83         0.1000   0.3610   0.26100   2.010
## 15    16    0.894  -2.17         0.9030   0.2820   0.19800   1.770
## 16    17    0.905  -2.10         0.9080   0.2240   0.15000   1.560
## 17    18    0.922  -2.04         0.9240   0.1810   0.11500   1.390
## 18    19    0.927  -2.01         0.9230   0.1480   0.08820   1.250
## 19    20    0.938  -1.93         0.9370   0.1220   0.06780   1.130
## 20    21    0.946  -1.86         0.9470   0.1020   0.05260   1.020
## 21    22    0.315  -3.12         0.2220   0.0856   0.04140   0.926
## 22    23    0.303  -3.96         0.1560   0.0728   0.03270   0.844
## 23    24    0.307  -2.92         0.1940   0.0623   0.02590   0.772
## 24    25    0.303  -2.84         0.1800   0.0538   0.02060   0.708
## 25    26    0.921  -1.63         0.8990   0.0467   0.01650   0.657
## 26    27    0.892  -1.61         0.8650   0.0408   0.01320   0.622
## 27    28    0.155  -1.93        -0.0202   0.0359   0.01040   0.591
## 28    29    0.199  -2.77        -0.0228   0.0317   0.00841   0.563
## 29    30    0.199  -2.69        -0.0203   0.0282   0.00684   0.537

plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")

0.7.2 vst

after vst normalisation

ps.varstab <- ps_vst(ps.inall, "Day")

data2 <- ps.varstab@otu_table@.Data %>% 
  as.data.frame()

rownames(data2) <- as.character(ps.varstab@sam_data$Description)
powers <-  c(seq(from = 1, to=10, by=0.5), seq(from = 11, to=20, by=1))
sft2 <-  pickSoftThreshold(data2, powerVector = powers, verbose = 5, networkType = "signed hybrid")

## pickSoftThreshold: will use block size 338.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 338 of 338
##    Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1    1.0    0.409 -1.05          0.704  50.100  45.60000  97.80
## 2    1.5    0.578 -1.07          0.895  30.900  27.00000  71.20
## 3    2.0    0.723 -1.15          0.942  20.300  17.90000  53.60
## 4    2.5    0.767 -1.24          0.871  14.000  12.00000  41.70
## 5    3.0    0.809 -1.31          0.909  10.000   7.86000  33.40
## 6    3.5    0.872 -1.31          0.958   7.380   5.39000  27.30
## 7    4.0    0.888 -1.30          0.968   5.580   3.76000  22.60
## 8    4.5    0.858 -1.36          0.910   4.320   2.67000  19.00
## 9    5.0    0.892 -1.36          0.943   3.400   1.94000  16.20
## 10   5.5    0.899 -1.38          0.954   2.720   1.41000  13.90
## 11   6.0    0.896 -1.39          0.962   2.210   1.07000  12.10
## 12   6.5    0.897 -1.33          0.962   1.820   0.84700  10.60
## 13   7.0    0.910 -1.30          0.970   1.520   0.65400   9.30
## 14   7.5    0.914 -1.31          0.966   1.280   0.51500   8.25
## 15   8.0    0.913 -1.28          0.954   1.090   0.41500   7.36
## 16   8.5    0.905 -1.28          0.937   0.938   0.34000   6.59
## 17   9.0    0.895 -1.30          0.913   0.813   0.27800   5.99
## 18   9.5    0.878 -1.33          0.900   0.709   0.23200   5.49
## 19  10.0    0.886 -1.34          0.894   0.624   0.19100   5.07
## 20  11.0    0.911 -1.31          0.928   0.491   0.12500   4.35
## 21  12.0    0.890 -1.29          0.894   0.396   0.08590   3.79
## 22  13.0    0.910 -1.26          0.918   0.325   0.06120   3.33
## 23  14.0    0.882 -1.25          0.856   0.272   0.04420   2.95
## 24  15.0    0.889 -1.19          0.864   0.230   0.03230   2.63
## 25  16.0    0.913 -1.18          0.893   0.198   0.02240   2.35
## 26  17.0    0.930 -1.16          0.912   0.172   0.01620   2.20
## 27  18.0    0.955 -1.19          0.947   0.151   0.01190   2.09
## 28  19.0    0.908 -1.21          0.896   0.134   0.00896   1.99
## 29  20.0    0.916 -1.21          0.906   0.120   0.00649   1.91

plot(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")

# WGCNA, as old and not supported
detachAllPackages()
library(WGCNA)


net3 <- WGCNA::blockwiseModules(data2,
                          power=5.5,
                          TOMType="signed",
                          networkType="signed hybrid",
                          nThreads=0)


mergedColors2 <-  WGCNA::labels2colors(net3$colors, colorSeq = c("salmon", "darkgreen", "cyan", "red", "blue", "plum"))

plotDendroAndColors(
  net3$dendrograms[[1]],
  mergedColors2[net3$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE,
  hang = 0.03,
  addGuide = TRUE,
  guideHang = 0.05)

unique(ps.inall@sam_data$Bag)

##  [1] D01 D03 D05 D07 D08 D10 D12 D13 D14 D15
## Levels: D01 D03 D05 D07 D08 D10 D12 D13 D14 D15

0.7.2.1 import libraries 2

library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(heatmaply)
library(WGCNA)
library(phyloseq)
library(ggtree)
library(tidyverse)
library(KneeArrower)

modules_of_interest = mergedColors2 %>% 
  unique()

module_df <- data.frame(
  asv = names(net3$colors),
  colors = mergedColors2
)

# module_df[module_df == "yellow"] <- "salmon"

submod <-  module_df %>%
  subset(colors %in% modules_of_interest)

row.names(module_df) = module_df$asv

subexpr = as.data.frame(t(data2))[submod$asv,]

submod_df <- data.frame(subexpr) %>%
  mutate(
    asv = row.names(.)
  ) %>%
  pivot_longer(-asv) %>%
  mutate(
    module = module_df[asv,]$colors
  )

submod_df <- submod_df %>% 
  mutate(name = gsub("\\_.*","",submod_df$name)) %>% 
  group_by(name, asv) %>% 
  summarise(value = mean(value), asv = asv, module = module) %>% 
  relocate(c(asv, name, value, module)) %>% 
  ungroup() %>% 
  mutate(module = as.factor(module))

p <- submod_df %>% 
  ggplot(., aes(x=name, y=value, group=asv)) +
  geom_line(aes(color = module),
            alpha = 0.2) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.position = "none") +
  facet_grid(rows = vars(module)) +
  labs(x = "day",
       y = "normalized abundance")

p + scale_color_manual(values = levels(submod_df$module)) +
  scale_x_discrete(labels=(c("D01"="3",
                                "D03"="14",
                                "D05"="28",
                                "D07"="49" ,
                                "D08"="63",
                                "D10"="91",
                                "D12"="119",
                                "D13"="140",
                                "D14"="161",
                                "D15"="182")
         ))

0.7.2.2 ordination

UNIFRAC - the sausage narrows down

ps.inall.col <- ps.inall
df <- module_df %>% 
  rename("id" = "asv")
df <- df %>% 
  dplyr::select(-"id") %>% 
  mutate(colors = as.factor(colors))
taxa <- as.data.frame(ps.inall@tax_table@.Data)
tx <- cbind(taxa, df)
tx$colors <- factor(tx$colors, levels = c("salmon", "darkgreen", "cyan", "red", "blue", "plum"))
tax_table(ps.inall.col) <- tax_table(as.matrix(tx))
ord <- ordinate(ps.inall.col, "NMDS", "unifrac")

## Run 0 stress 0.0685455 
## Run 1 stress 0.05568715 
## ... New best solution
## ... Procrustes: rmse 0.04076806  max resid 0.2160826 
## Run 2 stress 0.06784429 
## Run 3 stress 0.05628488 
## Run 4 stress 0.05568719 
## ... Procrustes: rmse 0.0001769095  max resid 0.0006062451 
## ... Similar to previous best
## Run 5 stress 0.07937976 
## Run 6 stress 0.07544596 
## Run 7 stress 0.05628489 
## Run 8 stress 0.05568719 
## ... Procrustes: rmse 0.00002559041  max resid 0.00008856774 
## ... Similar to previous best
## Run 9 stress 0.05568712 
## ... New best solution
## ... Procrustes: rmse 0.00002106782  max resid 0.00006934972 
## ... Similar to previous best
## Run 10 stress 0.06740137 
## Run 11 stress 0.06740136 
## Run 12 stress 0.06784421 
## Run 13 stress 0.06784436 
## Run 14 stress 0.06854546 
## Run 15 stress 0.07953205 
## Run 16 stress 0.05568719 
## ... Procrustes: rmse 0.00004255519  max resid 0.0001401031 
## ... Similar to previous best
## Run 17 stress 0.06740138 
## Run 18 stress 0.06728216 
## Run 19 stress 0.05568718 
## ... Procrustes: rmse 0.0000373923  max resid 0.0001312747 
## ... Similar to previous best
## Run 20 stress 0.07888487 
## *** Solution reached

plot_ordination(ps.inall.col, ord, type = "species", color = "colors") + 
  scale_color_manual(values =  c("salmon", "darkgreen", "cyan", "red", "blue", "plum")) +
  theme_bw() +
  theme(legend.position = "none")

Bray - the sausage even
Late stages on the left - with early clusters overlapping, late clusters seems separated What affects on the axis2? There is clearly some sort of pattern.

ord <- ordinate(ps.inall.col, "NMDS", "bray")

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.09036937 
## Run 1 stress 0.0903694 
## ... Procrustes: rmse 0.00001927814  max resid 0.00007086537 
## ... Similar to previous best
## Run 2 stress 0.09024752 
## ... New best solution
## ... Procrustes: rmse 0.003386048  max resid 0.01505777 
## Run 3 stress 0.09036938 
## ... Procrustes: rmse 0.003390695  max resid 0.01505826 
## Run 4 stress 0.09036937 
## ... Procrustes: rmse 0.003386035  max resid 0.01507802 
## Run 5 stress 0.09036937 
## ... Procrustes: rmse 0.003386432  max resid 0.01505862 
## Run 6 stress 0.1257473 
## Run 7 stress 0.09036937 
## ... Procrustes: rmse 0.003385691  max resid 0.01507976 
## Run 8 stress 0.1268919 
## Run 9 stress 0.09036944 
## ... Procrustes: rmse 0.003410952  max resid 0.01523599 
## Run 10 stress 0.09024752 
## ... New best solution
## ... Procrustes: rmse 0.00003844186  max resid 0.0001638862 
## ... Similar to previous best
## Run 11 stress 0.09036937 
## ... Procrustes: rmse 0.003380784  max resid 0.0150593 
## Run 12 stress 0.09036938 
## ... Procrustes: rmse 0.003379871  max resid 0.01508219 
## Run 13 stress 0.09036939 
## ... Procrustes: rmse 0.003385498  max resid 0.01512324 
## Run 14 stress 0.1358024 
## Run 15 stress 0.09024754 
## ... Procrustes: rmse 0.00003794727  max resid 0.0001665169 
## ... Similar to previous best
## Run 16 stress 0.1258672 
## Run 17 stress 0.1358024 
## Run 18 stress 0.09036937 
## ... Procrustes: rmse 0.003376715  max resid 0.01503468 
## Run 19 stress 0.126985 
## Run 20 stress 0.1360016 
## *** Solution reached

plot_ordination(ps.inall.col, ord, type = "species", color = "colors") + 
  scale_color_manual(values =  c("salmon", "darkgreen", "cyan", "red", "blue", "plum")) +
  theme_bw() +
  theme(legend.position = "none")

0.7.2.3 visualisation

Next are the same pictures for all the groups.

l_vst <- color_filt(ps.inall, submod_df)
l_vst

$blue \(blue\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 86 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 86 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 86 tips and 85 internal nodes ] refseq() DNAStringSet: [ 86 reference sequences ]

\(blue\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 86 89493 0 11387 1819 Avg#Reads 2556.94

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 86(100%) 86(100%) 86(100%) 85(98.84%) 79(91.86%) 50(58.14%) 7(8.14%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(blue\)heat \(blue\)heat_rel \(blue\)tree \(blue\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq314	Bacteria	Cyanobacteria	Vampirivibrionia	Obscuribacterales	Obscuribacteraceae	NA	NA
Seq23	Bacteria	Cyanobacteria	Vampirivibrionia	Obscuribacterales	Obscuribacteraceae	NA	NA
Seq65	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Sphingomonas	NA
Seq188	Bacteria	Proteobacteria	Alphaproteobacteria	Reyranellales	Reyranellaceae	Reyranella	soli
Seq334	Bacteria	Proteobacteria	Alphaproteobacteria	Reyranellales	Reyranellaceae	Reyranella	NA
Seq75	Bacteria	Proteobacteria	Alphaproteobacteria	Reyranellales	Reyranellaceae	Reyranella	NA
Seq167	Bacteria	Proteobacteria	Alphaproteobacteria	Reyranellales	Reyranellaceae	Reyranella	massiliensis
Seq273	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Kaistiaceae	Kaistia	NA
Seq131	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Pseudolabrys	NA
Seq138	Bacteria	Proteobacteria	Alphaproteobacteria	Micropepsales	Micropepsaceae	NA	NA
Seq104	Bacteria	Proteobacteria	Alphaproteobacteria	Micropepsales	Micropepsaceae	NA	NA
Seq129	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Altererythrobacter	NA
Seq28	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Devosiaceae	Devosia	NA
Seq360	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Mesorhizobium	NA
Seq292	Bacteria	Proteobacteria	Alphaproteobacteria	Caulobacterales	Caulobacteraceae	NA	NA
Seq166	Bacteria	Proteobacteria	Alphaproteobacteria	Caulobacterales	Hyphomonadaceae	Hirschia	NA
Seq39	Bacteria	Gemmatimonadota	Gemmatimonadetes	Gemmatimonadales	Gemmatimonadaceae	NA	NA
Seq30	Bacteria	Actinobacteriota	Thermoleophilia	Solirubrobacterales	Solirubrobacteraceae	Conexibacter	NA
Seq151	Bacteria	Actinobacteriota	Thermoleophilia	Solirubrobacterales	67-14	NA	NA
Seq250	Bacteria	Actinobacteriota	Thermoleophilia	Solirubrobacterales	67-14	NA	NA
Seq135	Bacteria	Myxococcota	Polyangia	Polyangiales	BIrii41	NA	NA
Seq46	Bacteria	Myxococcota	Polyangia	Polyangiales	BIrii41	NA	NA
Seq35	Bacteria	Myxococcota	Polyangia	Polyangiales	BIrii41	NA	NA
Seq106	Bacteria	Myxococcota	Polyangia	Polyangiales	Sandaracinaceae	NA	NA
Seq345	Bacteria	Bdellovibrionota	Oligoflexia	0319-6G20	NA	NA	NA
Seq198	Bacteria	Bdellovibrionota	Bdellovibrionia	Bdellovibrionales	Bdellovibrionaceae	Bdellovibrio	NA
Seq97	Bacteria	Bdellovibrionota	Bdellovibrionia	Bdellovibrionales	Bdellovibrionaceae	Bdellovibrio	NA
Seq183	Bacteria	Dependentiae	Babeliae	Babeliales	UBA12409	NA	NA
Seq424	Bacteria	Actinobacteriota	Actinobacteria	Propionibacteriales	Nocardioidaceae	Kribbella	NA
Seq136	Bacteria	Actinobacteriota	Acidimicrobiia	Microtrichales	Iamiaceae	Iamia	NA
Seq174	Bacteria	Actinobacteriota	Acidimicrobiia	Microtrichales	Iamiaceae	Iamia	NA
Seq43	Bacteria	Actinobacteriota	Actinobacteria	Micrococcales	Microbacteriaceae	Galbitalea	NA
Seq420	Bacteria	Actinobacteriota	Actinobacteria	Micromonosporales	Micromonosporaceae	Luedemannella	NA
Seq127	Bacteria	Actinobacteriota	Actinobacteria	Micromonosporales	Micromonosporaceae	Dactylosporangium	NA
Seq288	Bacteria	Actinobacteriota	Actinobacteria	Streptomycetales	Streptomycetaceae	NA	NA
Seq449	Bacteria	Actinobacteriota	Actinobacteria	Propionibacteriales	Propionibacteriaceae	Jiangella	NA
Seq91	Bacteria	Bacteroidota	Bacteroidia	Bacteroidales	NA	NA	NA
Seq149	Bacteria	Bacteroidota	Bacteroidia	Bacteroidetes VC2.1 Bac22	NA	NA	NA
Seq169	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	Sphingobacteriaceae	Mucilaginibacter	calamicampi
Seq81	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	NS11-12 marine group	NA	NA
Seq211	Bacteria	Bacteroidota	Bacteroidia	NA	NA	NA	NA
Seq119	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	env.OPS 17	NA	NA
Seq77	Bacteria	Spirochaetota	Spirochaetia	Spirochaetales	Spirochaetaceae	Salinispira	NA
Seq175	Bacteria	Spirochaetota	Spirochaetia	Spirochaetales	Spirochaetaceae	Spirochaeta 2	NA
Seq305	Bacteria	Chloroflexi	Anaerolineae	Anaerolineales	Anaerolineaceae	NA	NA
Seq379	Bacteria	Chloroflexi	Chloroflexia	Chloroflexales	Roseiflexaceae	NA	NA
Seq165	Bacteria	Chloroflexi	Chloroflexia	Chloroflexales	Roseiflexaceae	NA	NA
Seq327	Bacteria	Chloroflexi	Chloroflexia	Thermomicrobiales	JG30-KF-CM45	NA	NA
Seq252	Bacteria	Armatimonadota	Fimbriimonadia	Fimbriimonadales	Fimbriimonadaceae	NA	NA
Seq382	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Pedosphaerales	Pedosphaeraceae	NA	NA
Seq123	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Opitutales	Opitutaceae	Lacunisphaera	limnophila
Seq124	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Chthoniobacterales	Terrimicrobiaceae	Terrimicrobium	NA
Seq325	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Chthoniobacterales	Chthoniobacteraceae	LD29	NA
Seq272	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Verrucomicrobiales	Verrucomicrobiaceae	Roseimicrobium	gellanilyticum
Seq402	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Verrucomicrobiales	Rubritaleaceae	Luteolibacter	NA
Seq322	Bacteria	Planctomycetota	Planctomycetes	Pirellulales	Pirellulaceae	NA	NA
Seq356	Bacteria	Planctomycetota	Planctomycetes	Pirellulales	Pirellulaceae	NA	NA
Seq259	Bacteria	Planctomycetota	Planctomycetes	Pirellulales	Pirellulaceae	Pir4 lineage	NA
Seq279	Bacteria	Planctomycetota	Planctomycetes	Pirellulales	Pirellulaceae	Pir4 lineage	NA
Seq214	Bacteria	Planctomycetota	Planctomycetes	Planctomycetales	Rubinisphaeraceae	SH-PL14	NA
Seq209	Bacteria	Planctomycetota	Planctomycetes	Planctomycetales	Schlesneriaceae	Schlesneria	NA
Seq229	Bacteria	Acidobacteriota	Vicinamibacteria	Vicinamibacterales	Vicinamibacteraceae	NA	NA
Seq275	Bacteria	Acidobacteriota	Vicinamibacteria	Vicinamibacterales	Vicinamibacteraceae	NA	NA
Seq247	Bacteria	Acidobacteriota	Vicinamibacteria	Vicinamibacterales	NA	NA	NA
Seq316	Bacteria	Bacteroidota	Kapabacteria	Kapabacteriales	NA	NA	NA
Seq371	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq12	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq153	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq326	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq10	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	NA	NA
Seq82	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Terrimonas	NA
Seq462	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Pseudoflavitalea	NA
Seq150	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Flavitalea	NA
Seq290	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Edaphobaculum	NA
Seq339	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Edaphobaculum	NA
Seq120	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Taibaiella	NA
Seq233	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	211ds20	NA	NA
Seq133	Bacteria	Proteobacteria	Gammaproteobacteria	Steroidobacterales	Steroidobacteraceae	Steroidobacter	flavus
Seq349	Bacteria	Proteobacteria	Gammaproteobacteria	Steroidobacterales	Steroidobacteraceae	Steroidobacter	NA
Seq315	Bacteria	Proteobacteria	Gammaproteobacteria	Gammaproteobacteria Incertae Sedis	Unknown Family	Acidibacter	NA
Seq57	Bacteria	Proteobacteria	Gammaproteobacteria	Gammaproteobacteria Incertae Sedis	Unknown Family	Acidibacter	NA
Seq96	Bacteria	Proteobacteria	Gammaproteobacteria	R7C24	NA	NA	NA
Seq170	Bacteria	Proteobacteria	Gammaproteobacteria	Gammaproteobacteria Incertae Sedis	Unknown Family	Acidibacter	NA
Seq94	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Rhodanobacteraceae	Dokdonella	ginsengisoli
Seq114	Bacteria	Proteobacteria	Gammaproteobacteria	Diplorickettsiales	Diplorickettsiaceae	NA	NA
Seq235	Bacteria	Proteobacteria	Gammaproteobacteria	Legionellales	Legionellaceae	Legionella	NA

$cyan \(cyan\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 30 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 30 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 30 tips and 29 internal nodes ] refseq() DNAStringSet: [ 30 reference sequences ]

\(cyan\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 30 176073 336 15108 4911 Avg#Reads 5030.66

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 30(100%) 30(100%) 30(100%) 30(100%) 30(100%) 29(96.67%) 2(6.67%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(cyan\)heat \(cyan\)heat_rel \(cyan\)tree \(cyan\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq98	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	Burkholderia	telluris
Seq64	Bacteria	Proteobacteria	Alphaproteobacteria	Azospirillales	Inquilinaceae	Inquilinus	ginsengisoli
Seq4	Bacteria	Proteobacteria	Alphaproteobacteria	Azospirillales	Inquilinaceae	Inquilinus	NA
Seq244	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Pseudorhodoplanes	NA
Seq11	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Bradyrhizobium	NA
Seq74	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Bradyrhizobium	NA
Seq22	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Starkeya	NA
Seq45	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Allorhizobium	NA
Seq9	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Allorhizobium	NA
Seq29	Bacteria	Myxococcota	Polyangia	Nannocystales	Nannocystaceae	Nannocystis	NA
Seq239	Bacteria	Myxococcota	Polyangia	Polyangiales	Polyangiaceae	Labilithrix	NA
Seq50	Bacteria	Actinobacteriota	Actinobacteria	Corynebacteriales	Mycobacteriaceae	Mycobacterium	NA
Seq34	Bacteria	Actinobacteriota	Actinobacteria	Streptosporangiales	Streptosporangiaceae	Herbidospora	NA
Seq56	Bacteria	Firmicutes	Bacilli	Paenibacillales	Paenibacillaceae	Paenibacillus	NA
Seq25	Bacteria	Firmicutes	Bacilli	Bacillales	Bacillaceae	Terribacillus	NA
Seq3	Bacteria	Firmicutes	Bacilli	Bacillales	Bacillaceae	Bacillus	NA
Seq5	Bacteria	Firmicutes	Bacilli	Bacillales	Bacillaceae	Bacillus	NA
Seq13	Bacteria	Firmicutes	Bacilli	Bacillales	Planococcaceae	NA	NA
Seq7	Bacteria	Firmicutes	Bacilli	Bacillales	Planococcaceae	Solibacillus	NA
Seq19	Bacteria	Firmicutes	Bacilli	Bacillales	Planococcaceae	Solibacillus	NA
Seq17	Bacteria	Planctomycetota	Planctomycetes	Isosphaerales	Isosphaeraceae	Singulisphaera	NA
Seq324	Bacteria	Acidobacteriota	Acidobacteriae	Acidobacteriales	Acidobacteriaceae (Subgroup 1)	Edaphobacter	NA
Seq16	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq6	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq1	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq109	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq49	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Taibaiella	NA
Seq73	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Spirosomaceae	Dyadobacter	NA
Seq15	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Rhodanobacteraceae	Luteibacter	NA
Seq26	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Xanthomonadaceae	Luteimonas	NA

$darkgreen \(darkgreen\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 44 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 44 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 44 tips and 43 internal nodes ] refseq() DNAStringSet: [ 44 reference sequences ]

\(darkgreen\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 44 12360 0 2468 65 Avg#Reads 353.14

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 44(100%) 44(100%) 44(100%) 43(97.73%) 40(90.91%) 24(54.55%) 3(6.82%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(darkgreen\)heat \(darkgreen\)heat_rel \(darkgreen\)tree \(darkgreen\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq342	Archaea	Crenarchaeota	Nitrososphaeria	Nitrososphaerales	Nitrososphaeraceae	NA	NA
Seq196	Archaea	Crenarchaeota	Nitrososphaeria	Nitrososphaerales	Nitrososphaeraceae	Candidatus Nitrocosmicus	NA
Seq276	Archaea	Crenarchaeota	Nitrososphaeria	Nitrososphaerales	Nitrososphaeraceae	NA	NA
Seq419	Archaea	Crenarchaeota	Nitrososphaeria	Nitrososphaerales	Nitrososphaeraceae	NA	NA
Seq454	Bacteria	Proteobacteria	Gammaproteobacteria	Diplorickettsiales	Diplorickettsiaceae	Aquicella	NA
Seq329	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Sphingomonas	jaspsi
Seq332	Bacteria	Proteobacteria	Alphaproteobacteria	Rhodospirillales	Magnetospiraceae	NA	NA
Seq429	Bacteria	Proteobacteria	Alphaproteobacteria	Acetobacterales	Acetobacteraceae	NA	NA
Seq195	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Pseudorhodoplanes	NA
Seq144	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Pseudolabrys	NA
Seq369	Bacteria	Proteobacteria	Alphaproteobacteria	Micropepsales	Micropepsaceae	NA	NA
Seq362	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Altererythrobacter	NA
Seq428	Bacteria	Proteobacteria	Alphaproteobacteria	Caulobacterales	Caulobacteraceae	Phenylobacterium	NA
Seq347	Bacteria	Nitrospirota	Nitrospiria	Nitrospirales	Nitrospiraceae	Nitrospira	japonica
Seq328	Bacteria	Myxococcota	Polyangia	Haliangiales	Haliangiaceae	Haliangium	NA
Seq348	Bacteria	Myxococcota	Polyangia	Polyangiales	BIrii41	NA	NA
Seq308	Bacteria	Myxococcota	Polyangia	Polyangiales	BIrii41	NA	NA
Seq256	Bacteria	Myxococcota	Polyangia	Polyangiales	Polyangiaceae	Minicystis	NA
Seq302	Bacteria	Bdellovibrionota	Oligoflexia	0319-6G20	NA	NA	NA
Seq435	Bacteria	Actinobacteriota	Actinobacteria	Propionibacteriales	Propionibacteriaceae	Microlunatus	NA
Seq291	Bacteria	Actinobacteriota	Acidimicrobiia	IMCC26256	NA	NA	NA
Seq392	Bacteria	Actinobacteriota	Acidimicrobiia	Microtrichales	Ilumatobacteraceae	NA	NA
Seq321	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	NS11-12 marine group	NA	NA
Seq186	Bacteria	Firmicutes	Bacilli	Paenibacillales	Paenibacillaceae	Paenibacillus	NA
Seq204	Bacteria	Firmicutes	Bacilli	Bacillales	Planococcaceae	Paenisporosarcina	NA
Seq146	Bacteria	Spirochaetota	Spirochaetia	Spirochaetales	Spirochaetaceae	Spirochaeta 2	NA
Seq359	Bacteria	Patescibacteria	Saccharimonadia	Saccharimonadales	LWQ8	NA	NA
Seq269	Bacteria	Chloroflexi	Chloroflexia	Chloroflexales	Roseiflexaceae	NA	NA
Seq267	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Opitutales	Opitutaceae	Lacunisphaera	NA
Seq320	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Opitutales	Opitutaceae	Opitutus	NA
Seq203	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Chthoniobacterales	Terrimicrobiaceae	Terrimicrobium	NA
Seq497	Bacteria	Planctomycetota	Planctomycetes	Pirellulales	Pirellulaceae	NA	NA
Seq384	Bacteria	Planctomycetota	Planctomycetes	Planctomycetales	NA	NA	NA
Seq155	Bacteria	Acidobacteriota	Blastocatellia	Blastocatellales	Blastocatellaceae	NA	NA
Seq350	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	env.OPS 17	NA	NA
Seq224	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Terrimonas	NA
Seq409	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Edaphobaculum	NA
Seq337	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Taibaiella	NA
Seq494	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Amoebophilaceae	Candidatus Amoebophilus	NA
Seq277	Bacteria	Proteobacteria	Gammaproteobacteria	NA	NA	NA	NA
Seq191	Bacteria	Proteobacteria	Gammaproteobacteria	Gammaproteobacteria Incertae Sedis	Unknown Family	Acidibacter	NA
Seq335	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Xanthomonadaceae	Luteimonas	vadosa
Seq157	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Rhodanobacteraceae	Dokdonella	NA
Seq385	Bacteria	Proteobacteria	Gammaproteobacteria	Salinisphaerales	Solimonadaceae	NA	NA

$plum \(plum\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 30 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 30 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 30 tips and 29 internal nodes ] refseq() DNAStringSet: [ 30 reference sequences ]

\(plum\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 30 17172 0 3034 103 Avg#Reads 490.63

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 30(100%) 30(100%) 30(100%) 30(100%) 27(90%) 16(53.33%) 1(3.33%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(plum\)heat \(plum\)heat_rel \(plum\)tree \(plum\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq161	Bacteria	Proteobacteria	Gammaproteobacteria	Diplorickettsiales	Diplorickettsiaceae	Aquicella	NA
Seq181	Bacteria	Proteobacteria	Gammaproteobacteria	CCD24	NA	NA	NA
Seq207	Bacteria	Cyanobacteria	Vampirivibrionia	Obscuribacterales	Obscuribacteraceae	NA	NA
Seq219	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Sphingomonas	naasensis
Seq216	Bacteria	Proteobacteria	Alphaproteobacteria	Dongiales	Dongiaceae	Dongia	NA
Seq220	Bacteria	Proteobacteria	Alphaproteobacteria	Reyranellales	Reyranellaceae	Reyranella	NA
Seq90	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Pseudolabrys	NA
Seq271	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	NA	NA
Seq111	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Sphingobium	NA
Seq33	Bacteria	Proteobacteria	Alphaproteobacteria	Caulobacterales	Caulobacteraceae	Caulobacter	NA
Seq189	Bacteria	Actinobacteriota	Thermoleophilia	Solirubrobacterales	Solirubrobacteraceae	Solirubrobacter	NA
Seq85	Bacteria	Myxococcota	Polyangia	Polyangiales	BIrii41	NA	NA
Seq208	Bacteria	Actinobacteriota	Actinobacteria	Micrococcales	Microbacteriaceae	Galbitalea	NA
Seq218	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	Sphingobacteriaceae	Mucilaginibacter	NA
Seq215	Bacteria	Firmicutes	Bacilli	Bacillales	Bacillaceae	Bacillus	NA
Seq177	Bacteria	Fibrobacterota	Fibrobacteria	Fibrobacterales	Fibrobacteraceae	NA	NA
Seq375	Bacteria	Chloroflexi	Chloroflexia	Thermomicrobiales	JG30-KF-CM45	NA	NA
Seq430	Bacteria	Planctomycetota	Phycisphaerae	Phycisphaerales	Phycisphaeraceae	NA	NA
Seq368	Bacteria	Planctomycetota	Planctomycetes	Pirellulales	Pirellulaceae	NA	NA
Seq152	Bacteria	Planctomycetota	Planctomycetes	Gemmatales	Gemmataceae	Gemmata	NA
Seq304	Bacteria	Planctomycetota	Planctomycetes	Planctomycetales	NA	NA	NA
Seq141	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Terrimonas	NA
Seq343	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	NA	NA
Seq412	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Edaphobaculum	NA
Seq457	Bacteria	Proteobacteria	Gammaproteobacteria	Coxiellales	Coxiellaceae	Coxiella	NA
Seq130	Bacteria	Proteobacteria	Gammaproteobacteria	R7C24	NA	NA	NA
Seq346	Bacteria	Proteobacteria	Gammaproteobacteria	Salinisphaerales	Solimonadaceae	Alkanibacter	NA
Seq474	Bacteria	Proteobacteria	Gammaproteobacteria	Diplorickettsiales	Diplorickettsiaceae	NA	NA
Seq472	Bacteria	Proteobacteria	Gammaproteobacteria	Diplorickettsiales	Diplorickettsiaceae	NA	NA
Seq227	Bacteria	Proteobacteria	Gammaproteobacteria	Diplorickettsiales	Diplorickettsiaceae	NA	NA

$red \(red\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 75 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 75 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 75 tips and 74 internal nodes ] refseq() DNAStringSet: [ 75 reference sequences ]

\(red\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 75 66650 241 5339 1528 Avg#Reads 1904.29

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 75(100%) 75(100%) 75(100%) 74(98.67%) 74(98.67%) 66(88%) 11(14.67%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(red\)heat \(red\)heat_rel \(red\)tree \(red\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq117	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	Burkholderia	NA
Seq60	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Alcaligenaceae	NA	NA
Seq121	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Comamonadaceae	Variovorax	NA
Seq383	Bacteria	Proteobacteria	Alphaproteobacteria	Ferrovibrionales	Ferrovibrionaceae	Ferrovibrio	soli
Seq228	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	NA	NA
Seq178	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Beijerinckiaceae	Methylobacterium-Methylorubrum	NA
Seq171	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Beijerinckiaceae	Bosea	thiooxidans
Seq341	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Beijerinckiaceae	Bosea	NA
Seq89	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Tardiphaga	robiniae
Seq31	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Xanthobacteraceae	Starkeya	NA
Seq160	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Hyphomicrobiaceae	Hyphomicrobium	NA
Seq251	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Sphingopyxis	NA
Seq373	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Altererythrobacter	NA
Seq108	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Mesorhizobium	NA
Seq71	Bacteria	Actinobacteriota	Thermoleophilia	Solirubrobacterales	Solirubrobacteraceae	Conexibacter	NA
Seq338	Bacteria	Myxococcota	Polyangia	Polyangiales	Polyangiaceae	Pajaroellobacter	NA
Seq230	Bacteria	Myxococcota	Polyangia	Polyangiales	Polyangiaceae	Pajaroellobacter	NA
Seq88	Bacteria	Myxococcota	Polyangia	Polyangiales	Polyangiaceae	Sorangium	NA
Seq303	Bacteria	Actinobacteriota	Acidimicrobiia	Microtrichales	Ilumatobacteraceae	NA	NA
Seq93	Bacteria	Actinobacteriota	Actinobacteria	Micrococcales	Promicromonosporaceae	Promicromonospora	NA
Seq366	Bacteria	Actinobacteriota	Actinobacteria	Micrococcales	Microbacteriaceae	Leifsonia	NA
Seq134	Bacteria	Actinobacteriota	Actinobacteria	Micrococcales	Microbacteriaceae	NA	aoyamense
Seq411	Bacteria	Actinobacteriota	Actinobacteria	Micrococcales	Microbacteriaceae	Agromyces	NA
Seq115	Bacteria	Actinobacteriota	Actinobacteria	Corynebacteriales	Mycobacteriaceae	Mycobacterium	NA
Seq128	Bacteria	Actinobacteriota	Actinobacteria	Streptosporangiales	Streptosporangiaceae	Herbidospora	mongoliensis
Seq398	Bacteria	Actinobacteriota	Actinobacteria	Streptosporangiales	Streptosporangiaceae	Nonomuraea	NA
Seq78	Bacteria	Actinobacteriota	Actinobacteria	Streptosporangiales	Thermomonosporaceae	Actinocorallia	NA
Seq193	Bacteria	Actinobacteriota	Actinobacteria	Streptomycetales	Streptomycetaceae	Streptomyces	NA
Seq148	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	Sphingobacteriaceae	Pedobacter	NA
Seq307	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	NS11-12 marine group	NA	NA
Seq431	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	NS11-12 marine group	NA	NA
Seq205	Bacteria	Firmicutes	Bacilli	Paenibacillales	Paenibacillaceae	Paenibacillus	lautus
Seq86	Bacteria	Firmicutes	Bacilli	Paenibacillales	Paenibacillaceae	Paenibacillus	NA
Seq51	Bacteria	Firmicutes	Bacilli	Paenibacillales	Paenibacillaceae	Paenibacillus	NA
Seq18	Bacteria	Firmicutes	Bacilli	Paenibacillales	Paenibacillaceae	Paenibacillus	NA
Seq192	Bacteria	Firmicutes	Bacilli	Paenibacillales	Paenibacillaceae	Paenibacillus	NA
Seq266	Bacteria	Firmicutes	Bacilli	Paenibacillales	Paenibacillaceae	Paenibacillus	NA
Seq278	Bacteria	Firmicutes	Bacilli	Bacillales	Planococcaceae	Domibacillus	NA
Seq99	Bacteria	Firmicutes	Bacilli	Bacillales	Bacillaceae	Bacillus	NA
Seq53	Bacteria	Firmicutes	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	NA
Seq118	Bacteria	Firmicutes	Bacilli	Bacillales	Planococcaceae	Lysinibacillus	NA
Seq76	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Chthoniobacterales	Terrimicrobiaceae	Terrimicrobium	NA
Seq232	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Chthoniobacterales	Terrimicrobiaceae	Terrimicrobium	NA
Seq168	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Chthoniobacterales	Terrimicrobiaceae	Terrimicrobium	NA
Seq387	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Chthoniobacterales	Chthoniobacteraceae	Chthoniobacter	NA
Seq285	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Chthoniobacterales	Chthoniobacteraceae	Chthoniobacter	NA
Seq217	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Verrucomicrobiales	Verrucomicrobiaceae	Verrucomicrobium	spinosum
Seq190	Bacteria	Planctomycetota	Planctomycetes	Gemmatales	Gemmataceae	Gemmata	NA
Seq249	Bacteria	Planctomycetota	Planctomycetes	Isosphaerales	Isosphaeraceae	Singulisphaera	NA
Seq254	Bacteria	Acidobacteriota	Acidobacteriae	Acidobacteriales	Acidobacteriaceae (Subgroup 1)	Edaphobacter	NA
Seq67	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq176	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq268	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Pseudoflavitalea	NA
Seq32	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Pseudoflavitalea	NA
Seq84	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Pseudoflavitalea	NA
Seq459	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Niastella	NA
Seq40	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Niastella	hibisci
Seq173	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Niastella	NA
Seq246	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Niastella	NA
Seq112	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Flavitalea	NA
Seq126	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	ginsengihumi
Seq102	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	arvensicola
Seq309	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	soli
Seq70	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq262	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq301	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Taibaiella	NA
Seq572	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	Moraxellaceae	NA	NA
Seq436	Bacteria	Proteobacteria	Gammaproteobacteria	Steroidobacterales	Steroidobacteraceae	Steroidobacter	NA
Seq261	Bacteria	Proteobacteria	Gammaproteobacteria	Gammaproteobacteria Incertae Sedis	Unknown Family	Acidibacter	NA
Seq298	Bacteria	Proteobacteria	Gammaproteobacteria	Gammaproteobacteria Incertae Sedis	Unknown Family	Acidibacter	NA
Seq354	Bacteria	Proteobacteria	Gammaproteobacteria	Gammaproteobacteria Incertae Sedis	Unknown Family	Acidibacter	NA
Seq226	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Rhodanobacteraceae	Tahibacter	NA
Seq147	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Xanthomonadaceae	Xanthomonas	NA
Seq370	Bacteria	Proteobacteria	Gammaproteobacteria	Diplorickettsiales	Diplorickettsiaceae	NA	NA
Seq486	Bacteria	Proteobacteria	Gammaproteobacteria	NA	NA	NA	NA

$salmon \(salmon\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 73 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 73 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 73 tips and 72 internal nodes ] refseq() DNAStringSet: [ 73 reference sequences ]

\(salmon\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 73 115332 0 20534 779 Avg#Reads 3295.2

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 73(100%) 73(100%) 73(100%) 72(98.63%) 71(97.26%) 66(90.41%) 12(16.44%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(salmon\)heat \(salmon\)heat_rel \(salmon\)tree \(salmon\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq79	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Oxalobacteraceae	Massilia	armeniaca
Seq122	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Oxalobacteraceae	NA	NA
Seq179	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Oxalobacteraceae	Pseudoduganella	NA
Seq145	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Oxalobacteraceae	Pseudoduganella	eburnea
Seq154	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Oxalobacteraceae	Massilia	NA
Seq222	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Oxalobacteraceae	Massilia	NA
Seq159	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	Cupriavidus	NA
Seq8	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	Cupriavidus	NA
Seq14	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Alcaligenaceae	Achromobacter	NA
Seq54	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Comamonadaceae	NA	paradoxus
Seq95	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Comamonadaceae	Xylophilus	NA
Seq340	Bacteria	Proteobacteria	Gammaproteobacteria	Burkholderiales	Comamonadaceae	Rhizobacter	NA
Seq264	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Sphingomonas	mucosissima
Seq333	Bacteria	Proteobacteria	Alphaproteobacteria	Reyranellales	Reyranellaceae	Reyranella	NA
Seq158	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Beijerinckiaceae	Microvirga	NA
Seq248	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Allorhizobium	NA
Seq245	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Devosiaceae	Devosia	neptuniae
Seq185	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Shinella	NA
Seq156	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	NA	NA
Seq42	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	NA	NA
Seq286	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Ensifer	NA
Seq21	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Allorhizobium	NA
Seq107	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Neorhizobium	NA
Seq231	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	Allorhizobium	azooxidifex
Seq140	Bacteria	Proteobacteria	Alphaproteobacteria	Rhizobiales	Rhizobiaceae	NA	NA
Seq378	Bacteria	Proteobacteria	Alphaproteobacteria	Caulobacterales	Caulobacteraceae	Phenylobacterium	mobile
Seq374	Bacteria	Myxococcota	Polyangia	Haliangiales	Haliangiaceae	Haliangium	NA
Seq257	Bacteria	Myxococcota	Polyangia	Polyangiales	Polyangiaceae	Pajaroellobacter	NA
Seq223	Bacteria	Actinobacteriota	Actinobacteria	Micrococcales	Promicromonosporaceae	Cellulosimicrobium	NA
Seq83	Bacteria	Actinobacteriota	Actinobacteria	Micrococcales	Microbacteriaceae	Microbacterium	NA
Seq68	Bacteria	Actinobacteriota	Actinobacteria	Corynebacteriales	Mycobacteriaceae	Mycobacterium	NA
Seq293	Bacteria	Actinobacteriota	Actinobacteria	Glycomycetales	Glycomycetaceae	Glycomyces	NA
Seq72	Bacteria	Actinobacteriota	Actinobacteria	Streptomycetales	Streptomycetaceae	Streptomyces	NA
Seq101	Bacteria	Actinobacteriota	Actinobacteria	Streptomycetales	Streptomycetaceae	Streptomyces	NA
Seq265	Bacteria	Actinobacteriota	Actinobacteria	Streptomycetales	Streptomycetaceae	Streptomyces	NA
Seq260	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	NA	NA	NA
Seq180	Bacteria	Bacteroidota	Bacteroidia	Sphingobacteriales	Sphingobacteriaceae	Pedobacter	panaciterrae
Seq142	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Verrucomicrobiales	Verrucomicrobiaceae	Roseimicrobium	NA
Seq163	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq253	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq116	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq184	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Microscillaceae	Ohtaekwangia	NA
Seq52	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Flavobacteriaceae	Flavobacterium	NA
Seq187	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Weeksellaceae	Chryseobacterium	ginsenosidimutans
Seq172	Bacteria	Bacteroidota	Bacteroidia	Flavobacteriales	Weeksellaceae	Chryseobacterium	NA
Seq283	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Flavitalea	NA
Seq113	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Niastella	NA
Seq110	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Niastella	NA
Seq132	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Pseudoflavitalea	NA
Seq37	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq162	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq182	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	humicola
Seq92	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq80	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	pinensis
Seq2	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq221	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq20	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq282	Bacteria	Bacteroidota	Bacteroidia	Chitinophagales	Chitinophagaceae	Chitinophaga	NA
Seq38	Bacteria	Bacteroidota	Bacteroidia	Cytophagales	Spirosomaceae	Dyadobacter	NA
Seq237	Bacteria	Proteobacteria	Gammaproteobacteria	Steroidobacterales	Steroidobacteraceae	Steroidobacter	agariperforans
Seq234	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Xanthomonadaceae	Pseudoxanthomonas	NA
Seq59	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Xanthomonadaceae	Stenotrophomonas	NA
Seq105	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Xanthomonadaceae	Stenotrophomonas	NA
Seq344	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Xanthomonadaceae	Lysobacter	NA
Seq41	Bacteria	Proteobacteria	Gammaproteobacteria	Xanthomonadales	Xanthomonadaceae	Lysobacter	NA
Seq243	Bacteria	Proteobacteria	Gammaproteobacteria	NA	NA	NA	NA
Seq406	Bacteria	Proteobacteria	Gammaproteobacteria	Legionellales	Legionellaceae	Legionella	NA
Seq139	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	Pseudomonadaceae	Pseudomonas	NA
Seq24	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	Pseudomonadaceae	Pseudomonas	NA
Seq27	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	Pseudomonadaceae	Pseudomonas	NA
Seq87	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	Pseudomonadaceae	Pseudomonas	NA
Seq55	Bacteria	Proteobacteria	Gammaproteobacteria	Pseudomonadales	Pseudomonadaceae	Pseudomonas	NA
Seq201	Bacteria	Proteobacteria	Gammaproteobacteria	Enterobacterales	Enterobacteriaceae	Klebsiella	NA

0.7.3 mdp and other indices for clusters

I can’t do stat mdp support for clusters, but the early clusters(salmon, cyan) have lower mdp than the later ones.

external_empty_dataframe <- data.frame(cluster=factor(),
                                       mdp=numeric(),
                                       richness=numeric(),
                                       stringsAsFactors = FALSE)

for (i in names(l_vst)) {
  ps.in <- l_vst[[i]][["ps"]]
  m <- ps.in@otu_table@.Data %>% 
    colSums() %>% 
    as.data.frame() %>% 
    as.matrix() %>% 
    t()
  
  mdp_index <- picante::mpd(m, cophenetic(ps.in@phy_tree)) 
  ricness_index <- ps.in@tax_table %>% 
    rownames() %>% 
    length()
  d <- data.frame(cluster = i,
                  mdp = mdp_index,
                  richness = ricness_index)
  external_empty_dataframe <- rbind(external_empty_dataframe, d)
}

external_empty_dataframe %>% 
  dplyr::arrange(mdp)

##     cluster      mdp richness
## 1      cyan 1.567685       30
## 2    salmon 1.684152       73
## 3       red 1.748594       75
## 4      plum 1.779306       30
## 5      blue 1.798298       86
## 6 darkgreen 2.203539       44

colSums(l_vst$salmon$ps@otu_table@.Data)

##  Seq79 Seq122 Seq179 Seq145 Seq154 Seq222 Seq159   Seq8  Seq14  Seq54  Seq95 
##   1570    975    549    765    703    390    669   9841   6482   2589   1262 
## Seq340 Seq264 Seq333 Seq158 Seq248 Seq245 Seq185 Seq156  Seq42 Seq286  Seq21 
##    195    304    200    669    337    341    516    688   3429    261   5170 
## Seq107 Seq231 Seq140 Seq378 Seq374 Seq257 Seq223  Seq83  Seq68 Seq293  Seq72 
##   1064    363    798    158    160    310    390   1469   1941    249   1822 
## Seq101 Seq265 Seq260 Seq180 Seq142 Seq163 Seq253 Seq116 Seq184  Seq52 Seq187 
##   1179    301    309    542    791    638    317   1038    523   2663    512 
## Seq172 Seq283 Seq113 Seq110 Seq132  Seq37 Seq162 Seq182  Seq92  Seq80   Seq2 
##    572    266   1047   1054    875   3566    653    532   1311   1524  16855 
## Seq221  Seq20 Seq282  Seq38 Seq237 Seq234  Seq59 Seq105 Seq344  Seq41 Seq243 
##    392   5184    266   3539    350    357   2372   1105    189   3465    345 
## Seq406 Seq139  Seq24  Seq27  Seq87  Seq55 Seq201 
##    137    799   4452   4271   1373   2582    457

0.8 Add a metadata

list.files("meta/")

## [1] "cell_realtime_stat.xlsx" "cell_resp_ch_stat.xlsx" 
## [3] "period_legend.xlsx"

realtime.data <- readxl::read_excel("meta/cell_realtime_stat.xlsx")
period.data <- readxl::read_excel("meta/period_legend.xlsx")
resp.data <- readxl::read_excel("meta/cell_resp_ch_stat.xlsx")
realtime.data

## # A tibble: 36 × 29
##    id        day contr…¹ CELL_…² CELL_…³ CELL_…⁴ CELL_…⁵ CELL_…⁶ CELL_…⁷ CELL_…⁸
##    <chr>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 G01-1-1     3    15.6    33.2    31.1    33.6    36.6    33.8    38.0    33.0
##  2 G01-1-2     3    15.5    33.3    31.3    32.9    36.2    33.8    41      32.7
##  3 G01-1-3     3    15.5    33.6    31.2    32.9    35.8    33.3    38.2    32.9
##  4 G01-2-1     3    17.1    36.0    31.1    33.1    36.9    36.4    38.8    33.8
##  5 G01-2-2     3    17.3    35.6    31.0    32.8    37.8    35.7    40.1    33.6
##  6 G01-2-3     3    17.1    36      31.1    32.8    38.1    35.2    37.6    33.6
##  7 G01-3-1     3    16.4    35.4    31.6    34.9    39.1    35.3    40.9    33.2
##  8 G01-3-2     3    16.5    35.6    31.3    34.4    37.9    34.6    38.9    33.1
##  9 G01-3-3     3    16.4    35.8    31.5    34.4    41.6    35.0    36.5    33.4
## 10 G05-1-1    28    18.2    26.3    29.1    29.4    28.2    28.4    37.4    34.5
## # … with 26 more rows, 19 more variables: CELL_193122 <dbl>, CELL_73229 <dbl>,
## #   CELL_47814 <dbl>, CELL_163125 <dbl>, CELL_73266 <dbl>, CELL_88582 <dbl>,
## #   CELL_63583 <dbl>, CELL_14199 <dbl>, CELL_95616 <dbl>, CELL_63504 <dbl>,
## #   CELL_08643 <chr>, CELL_199599 <dbl>, CELL_01426 <dbl>, CELL_71601 <dbl>,
## #   CELL_45099 <dbl>, CELL_191900 <dbl>, CELL_99463 <dbl>, CELL_74579 <dbl>,
## #   CELL_183403 <dbl>, and abbreviated variable names ¹​contr_16S, ²​CELL_172283,
## #   ³​CELL_203163, ⁴​CELL_83325, ⁵​CELL_188413, ⁶​CELL_109631, ⁷​CELL_188119, …

0.8.1 Realtime data

impute.mean <- function(x) replace(x, is.na(x), mean(x, na.rm = TRUE))

realtime_data <- realtime.data %>% 
  mutate(CELL_08643 = as.numeric(CELL_08643)) %>% 
  group_by(day) %>%
  mutate(nice_cell = impute.mean(CELL_08643)) %>% 
  mutate(day = as.factor(day))

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

realtime_zjena <- readxl::read_excel("cellulases_gene_expression(1).xlsx")
geom_mean <-function(x){exp(mean(log(x)))}

realtime_zjena_geom <- realtime_zjena %>% 
  mutate(repeats = paste0(realtime_zjena$week, "-", rep(1:3, 3, each=3))) %>% 
  relocate(repeats, 1) %>% 
  group_by(repeats) %>% 
  summarise_if(is.numeric, geom_mean) %>% 
  mutate(week = as.factor(week)) %>% 
  arrange(week)

UPGMA cellulase(distance - bray).
I don’t really know what that means. Depending on the distance, you get very different clusters.

realtime_matrix <- realtime_zjena_geom %>% 
  column_to_rownames("repeats") %>% 
  select_if(is.numeric) %>% 
  as.matrix()

hcl <- hclust(vegan::vegdist(t(realtime_matrix), method="bray"), "average")
plot(hcl)

And here is the Euclidean distance.
Nothing to compare with the previous picture

hcl <- hclust(vegan::vegdist(t(realtime_matrix), method="euclidian"), "average")
plot(hcl)

Cellulases corplot.
Correlations only

m = cor(realtime_matrix)
corrplot::corrplot(m)

Only significant correlations . In general, there are clusters, but the correlations are not significant (with some exceptions)

cor_test_mat <- psych::corr.test(realtime_matrix)$p
corrplot::corrplot(m, p.mat = cor_test_mat, method = 'circle', type = 'lower', insig='blank',
         order = 'AOE', diag = FALSE)$corrPos -> p1

# text(p1$x, p1$y, round(p1$corr, 2))

The same, but the matrix is already logarithmased

cor_test_mat <- psych::corr.test(log(realtime_matrix))
corrplot::corrplot(cor_test_mat$r, p.mat = cor_test_mat$p, method = 'circle', type = 'lower', insig='blank',
         order = 'AOE', diag = FALSE)$corrPos -> p1

# text(p1$x, p1$y, round(p1$corr, 2))

add phylogenic tree (mafft - iqtree(ML))

realtime_tree <- ape::read.tree("al_chz.fasta.contree")
plot(realtime_tree)

library(tidyverse)

realtime_data %>% 
  select(c("day", "id", "contr_16S", "CELL_172283")) %>% 
  mutate(bio_repl = gsub("-[1-3]$", "", id)) %>% 
  group_by(bio_repl, day) %>% 
  summarise(contr_tmean = mean(contr_16S),
            data_tmean = mean(CELL_172283)) %>% 
  mutate(dCt = data_tmean - contr_tmean)

## `summarise()` has grouped output by 'bio_repl'. You can override using the
## `.groups` argument.

## # A tibble: 12 × 5
## # Groups:   bio_repl [12]
##    bio_repl day   contr_tmean data_tmean   dCt
##    <chr>    <fct>       <dbl>      <dbl> <dbl>
##  1 G01-1    3            15.5       33.4 17.8 
##  2 G01-2    3            17.2       35.9 18.7 
##  3 G01-3    3            16.4       35.6 19.2 
##  4 G05-1    28           18.4       26.2  7.81
##  5 G05-2    28           16.5       26.2  9.73
##  6 G05-3    28           18.5       32.0 13.6 
##  7 G10-1    91           17.4       26.2  8.72
##  8 G10-2    91           16.3       24.1  7.79
##  9 G10-3    91           17.5       26.3  8.81
## 10 G14-1    161          18.0       28.1 10.1 
## 11 G14-2    161          16.5       27.0 10.6 
## 12 G14-3    161          16.9       25.4  8.41

resp_data <- resp.data %>% 
  group_by(day) %>%
  mutate(control = impute.mean(control)) %>% 
  mutate(straw = impute.mean(straw)) 

resp_data

## # A tibble: 164 × 3
## # Groups:   day [27]
##      day control straw
##    <dbl>   <dbl> <dbl>
##  1     0    250   250 
##  2     3    285.  723.
##  3     3    263.  832.
##  4     3    131.  657.
##  5     3    307.  525.
##  6     3    241.  744.
##  7     3    245.  701.
##  8     3    245.  657.
##  9     7    274.  690.
## 10     7    296.  690.
## # … with 154 more rows

0.8.2 Knee plot

Soil resperation - median(straw)/median(control)

You can use a third-party package (KneeArrower) to understand where the knee_plot breaks. For the elimination of repetitions let’s take the median.   I don’t know how this package works (looks for a derivative, but I don’t know how it smoothes it out, it’s math), but it says that the breaking is more likely at 60-80 days.

(I corrected the errors - got closer to your data)

https://github.com/agentlans/KneeArrower - here’s where you can read about math

# resp_data %>% 
#   filter(!day == 0) %>% 
#   group_by(day) %>% 
#   summarise(median_control = median(control),
#             median_straw = median(straw)) %>% 
#   mutate(rel = median_straw/median_control) %>% 
#   ggplot() +
#   geom_point(aes(x = day, y = rel))

resp_median <- resp_data %>% 
  filter(!day == 0) %>% 
  group_by(day) %>% 
  summarise(median_control = median(control),
            median_straw = median(straw)) %>% 
  mutate(rel = median_straw/median_control)


thresholds <- c(0.25, 0.5, 0.75, 1)

# Find cutoff points at each threshold
cutoff.points <- lapply(thresholds, function(i) {
  findCutoff(resp_median$day, resp_median$rel, method="first", i)
})

x.coord <- sapply(cutoff.points, function(p) p$x)
y.coord <- sapply(cutoff.points, function(p) p$y)

# Plot the cutoff points on the scatterplot
plot(resp_median$day, resp_median$rel, pch=20, col="gray")
points(x.coord, y.coord, col="red", pch=20)
text(x.coord, y.coord, labels=thresholds, pos=4, col="red")

period_data <- period.data %>%
  mutate(bag_id = as.factor(bag_id))

period_data

## # A tibble: 15 × 2
##    bag_id   day
##    <fct>  <dbl>
##  1 1          3
##  2 2          7
##  3 3         14
##  4 4         21
##  5 5         28
##  6 6         35
##  7 7         49
##  8 8         63
##  9 9         77
## 10 10        91
## 11 11       105
## 12 12       119
## 13 13       140
## 14 14       161
## 15 15       182

0.8.3 alpha/respiration.

Well, I don’t really understand what to do next - attach the clusters to this picture?

resp_median_bags <- resp_median %>% 
  left_join(period.data, by="day") %>% 
  mutate(bag_id = as.factor(bag_id))

ps.f.r <- rarefy_even_depth(ps.f)

estimate_richness(ps.f.r) %>% 
  rownames_to_column("ID") %>% 
  mutate(bag_id = as.factor(
    as.numeric(
      gsub("\\..+","",
           gsub("straw\\.16s\\.D","", ID)
           )
      )
    )
    ) %>% 
  group_by(bag_id) %>% 
  summarise(Observed = mean(Observed),
            Shannon = mean(Shannon),
            InvSimpson = mean(InvSimpson)) %>% 
  left_join(period_data, by="bag_id") %>% 
  left_join(resp_median_bags, by="bag_id") %>% 
  mutate(
    Observed_scaled = scale(Observed),
    Shannon_scaled = scale(Shannon),
    InvSimpson_scaled = scale(InvSimpson),
    Respiration_scaled = scale(rel)
  ) %>% 
  select(c(bag_id, Observed_scaled, Shannon_scaled, InvSimpson_scaled, Respiration_scaled)) %>% 
  pivot_longer(c("Observed_scaled", "Shannon_scaled", "InvSimpson_scaled", "Respiration_scaled")) %>% 
  ggplot(aes(y = value, x = bag_id, group = name)) +
  geom_line(aes(color = name),
            alpha = 0.8) +
  theme_bw()

Correlation is negative, weak, not particularly reliable, and only Pearson (which works for a normal distribution (we have a time series, we need a spearman for good measure))

alpha_resp <- phyloseq::estimate_richness(ps.f.r) %>% 
  rownames_to_column("ID") %>% 
  mutate(bag_id = as.factor(
    as.numeric(
      gsub("\\..+","",
           gsub("straw\\.16s\\.D","", ID)
           )
      )
    )
    ) %>% 
  group_by(bag_id) %>% 
  summarise(Observed = mean(Observed),
            Shannon = mean(Shannon),
            InvSimpson = mean(InvSimpson)) %>% 
  left_join(period_data, by="bag_id") %>% 
  left_join(resp_median_bags, by="bag_id") %>% 
  mutate(
    Observed_scaled = scale(Observed),
    Shannon_scaled = scale(Shannon),
    InvSimpson_scaled = scale(InvSimpson),
    Respiration_scaled = scale(rel)
  ) %>% 
  select(c(bag_id, Observed_scaled, Shannon_scaled, InvSimpson_scaled, Respiration_scaled))

cor.test(alpha_resp$Observed_scaled, alpha_resp$Respiration_scaled, method = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  alpha_resp$Observed_scaled and alpha_resp$Respiration_scaled
## S = 230, p-value = 0.2629
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.3939394

cor.test(alpha_resp$Shannon_scaled, alpha_resp$Respiration_scaled, method = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  alpha_resp$Shannon_scaled and alpha_resp$Respiration_scaled
## S = 224, p-value = 0.3128
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.3575758

cor.test(alpha_resp$InvSimpson_scaled, alpha_resp$Respiration_scaled, method = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  alpha_resp$InvSimpson_scaled and alpha_resp$Respiration_scaled
## S = 226, p-value = 0.2956
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.369697

cor.test(alpha_resp$Observed_scaled, alpha_resp$Respiration_scaled, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  alpha_resp$Observed_scaled and alpha_resp$Respiration_scaled
## t = -2.7866, df = 8, p-value = 0.02368
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9234071 -0.1293521
## sample estimates:
##        cor 
## -0.7018198

cor.test(alpha_resp$Shannon_scaled, alpha_resp$Respiration_scaled, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  alpha_resp$Shannon_scaled and alpha_resp$Respiration_scaled
## t = -2.8622, df = 8, p-value = 0.02108
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9261458 -0.1479056
## sample estimates:
##        cor 
## -0.7112926

cor.test(alpha_resp$InvSimpson_scaled, alpha_resp$Respiration_scaled, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  alpha_resp$InvSimpson_scaled and alpha_resp$Respiration_scaled
## t = -2.2589, df = 8, p-value = 0.05381
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.900037185  0.009178016
## sample estimates:
##        cor 
## -0.6240545

0.9 Export all phyloseq object in one table

In the table there is a column cluster - “exl” in it means that this is the part of the dataset that did not go to WGCNA - majors in single samples(OMs).

ps.m <- phyloseq::psmelt(ps.f)

ps.m <- ps.m %>% 
  mutate_if(is.character, as.factor) 

ps.data.out <- ps.m %>%
  select(-Group) %>% 
  pivot_wider(names_from = c(Day, Description, Sample), values_from = Abundance, values_fill = 0) 

#create empty dataframe with columnnames
external_empty_dataframe <- data.frame(OTU=factor(), cluster=factor(), stringsAsFactors = FALSE)

for (i in names(l_vst)) {
  a <- taxa_names(l_vst[[i]][["ps"]])
  b <- rep(i, length(a))
  d <- data.frame(OTU = as.factor(a),
                  cluster = as.factor(b))
  external_empty_dataframe <- rbind(external_empty_dataframe, d)
}

clusters.otu.df <- external_empty_dataframe
# add exl taxa -- taxa "exl"
d <- data.frame(OTU = as.factor(ps.exl.taxa),
                cluster = as.factor(rep("exl", length(ps.exl.taxa))))

clusters.otu.df <- rbind(clusters.otu.df, d)
ps.data.out.exl <- left_join(clusters.otu.df, ps.data.out, by="OTU")
# write.table(ps.data.out.exl, file = "ps.data.out.tsv", sep = "\t")
ps.data.out.exl

0.9.1 Several pictures on the dynamics of representation of individual phyla

Bdellovibrionota - predators, an indicator of a developed community

ps.m %>% 
  filter(Phylum == "Bdellovibrionota") %>% 
  group_by(Description, Day) %>% 
  summarise(Bs = sum(Abundance)) %>% 
  ggplot() +
  geom_boxplot(aes(x = Day, y = Bs)) +
  theme_bw()

## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.

Myxococcota - sort of the same, but as we know can be cellulotic(facultative predators)

ps.m %>% 
  filter(Phylum == "Myxococcota") %>% 
  group_by(Description, Day) %>% 
  summarise(Bs = sum(Abundance)) %>% 
  ggplot() +
  geom_boxplot(aes(x = Day, y = Bs)) +
  theme_bw()

## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.

Археи появляются тоже на поздних стадиях - вообще я бы хотел бы опять развить тему важности азотного метаболизма на поздних стадиях разложения.
Оч хочется метагеном, но не этот.
Кроме того хочется отметить, что эти минорные группы возникают на D12

ps.m %>% 
  filter(Phylum == "Crenarchaeota") %>% 
  group_by(Description, Day) %>% 
  summarise(Bs = sum(Abundance)) %>% 
  ggplot() +
  geom_boxplot(aes(x = Day, y = Bs)) +
  theme_bw()

## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.

Gammaproteobacteria - they are a third of the cluster blue

ps.m %>% 
  filter(Class == "Gammaproteobacteria") %>% 
  group_by(Description, Day) %>% 
  summarise(Bs = sum(Abundance)) %>% 
  ggplot() +
  geom_boxplot(aes(x = Day, y = Bs)) +
  theme_bw()

## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.

0.9.2 Boxplots by phyla

All phyla - the representation is logarithmic on the basis of 2
Separate points are sums of absolute values of reads by days
The sums are logarithmic, but not the individual phylotypes

#select only major phylums
top_phylum <- ps.m %>% 
  count(Phylum) %>% 
  arrange(desc(n)) %>% 
  top_n(10) %>% 
  pull(Phylum)

## Selecting by n

ps.m %>% 
  filter(Phylum %in% top_phylum) %>% 
  mutate(
    Phylum = as.character(Phylum),
    Class = as.character(Class),
    phylum = ifelse(Phylum == "Proteobacteria", Class, Phylum)
  ) %>% 
  group_by(Description, Day, phylum) %>% 
  filter(!is.na(phylum)) %>% 
  summarise(Bs = log2(sum(Abundance))) %>% 
  ggplot(aes(x = Day, y = Bs)) +
  geom_boxplot(fill="#4DBBD5B2", alpha=0.4) +
  theme_bw() +
  facet_wrap(~ phylum)

## `summarise()` has grouped output by 'Description', 'Day'. You can override
## using the `.groups` argument.

## Warning: Removed 46 rows containing non-finite values (stat_boxplot).

0.9.3 resperation and phylotypes correlations

little_res <- resp_median_bags %>% 
  filter(day %in% ps.varstab@sam_data$Day)

ps.varstab.merged <- phyloseq::merge_samples(ps.varstab, "Day")

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'tidytree'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'tidytree'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'tidytree'

## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'

## Also defined by 'tidytree'

otus_cor <- ps.varstab.merged@otu_table %>% 
  as.data.frame()

s1 <- sapply(colnames(otus_cor), function(x) {
  cor <- cor.test(otus_cor[[x]], little_res$rel, method="spearman", exact = FALSE)
  c1 <- cor$p.value
  c2 <- cor$estimate %>% unname()
  return(list("p.value"=c1, "rho"=c2))})

s1 <- s1 %>% 
  t() %>% 
  as.data.frame() %>% 
  unnest(cols = c(p.value, rho), .id = "id") %>% 
  column_to_rownames("id") %>% 
  mutate(p.adj = p.adjust(p.value, method = "BH", n = length(p.value)))


taxa <- as.data.frame(ps.varstab@tax_table@.Data) %>% 
  rownames_to_column("ID")

s1 %>%
  rownames_to_column("ID") %>% 
  left_join(taxa) %>% 
  arrange(rho) %>%
  filter(Phylum %in% c("Proteobacteria", "Acidobacteriota", "Actinobacteriota","Chloroflexi", "Planctomycetota", "Bacteroidota", "Verrucomicrobiota", "Firmicutes")) %>% 
  mutate(Rank = ifelse(Phylum == "Proteobacteria", Class, Phylum)) %>%
  ggplot(aes(x=rho, color=Rank, fill=Rank)) +
  geom_histogram(alpha=0.3, bins = 30) +
  theme_minimal() +
  facet_wrap(~ Rank) +
  theme(legend.position="none")

## Joining, by = "ID"

s1.otus <- s1 %>% 
  rownames_to_column("OTU")

ps.data.out.exl %>% 
  select_if(Negate(is.integer)) %>% 
  select(-Bag) %>% 
  distinct() %>% 
  left_join(s1.otus, by="OTU") %>%
  filter(rho >= 0) %>%
  arrange(p.adj) %>% 
  head()

##     OTU cluster  Kingdom         Phylum               Class           Order
## 1 Seq24  salmon Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales
## 2 Seq27  salmon Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales
## 3  Seq2  salmon Bacteria   Bacteroidota         Bacteroidia Chitinophagales
## 4 Seq20  salmon Bacteria   Bacteroidota         Bacteroidia Chitinophagales
## 5 Seq59  salmon Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales
## 6 Seq41  salmon Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales
##             Family            Genus Species     p.value       rho     p.adj
## 1 Pseudomonadaceae      Pseudomonas    <NA> 0.005446342 0.8000947 0.1048480
## 2 Pseudomonadaceae      Pseudomonas    <NA> 0.005579938 0.7987953 0.1048480
## 3 Chitinophagaceae     Chitinophaga    <NA> 0.015800596 0.7333333 0.1372353
## 4 Chitinophagaceae     Chitinophaga    <NA> 0.023456474 0.7026122 0.1618018
## 5 Xanthomonadaceae Stenotrophomonas    <NA> 0.052195671 0.6273323 0.1745491
## 6 Xanthomonadaceae       Lysobacter    <NA> 0.066688000 0.6000000 0.1753438

0.10 ITS

same analysis again, but without

library(DESeq2)

ps_its <- readRDS("../../its/d2/ps_its")

sample.data <- ps_its@sam_data %>% 
  data.frame() %>% 
  mutate(Group = if_else(Day %in% c("D01", "D03", "D05"), "early", 
                         if_else(Day %in% c("D07", "D08","D10"), "middle", "late"))) %>% 
  mutate(Group = factor(Group, levels=c("early", "middle","late"))) %>% 
  mutate(Day = Day %>% 
           forcats::fct_recode( "3" = "D01",
                                "14" = "D03",
                                "28" = "D05",
                                "49" = "D07" ,
                                "63" = "D08",
                                "91" = "D10",
                                "119" = "D12",
                                "140" = "D13",
                                "161" = "D14",
                                "182" = "D15")
         ) %>% 
  phyloseq::sample_data()

sample_data(ps_its) <- sample.data

amp_its <- phyloseq_to_ampvis2(ps_its)

ps.its.inall <- phyloseq::filter_taxa(ps_its, function(x) sum(x > 10) > (0.1*length(x)), TRUE)
ps.its.inall <- prune_taxa(taxa_sums(ps.its.inall) > 0, ps.its.inall)
Biostrings::writeXStringSet(ps.its.inall@refseq, file="ref_inall.fasta")
Biostrings::writeXStringSet(ps_its@refseq, file="ref.fasta")

tree <- ape::read.tree("../../its/d2/al.fasta.contree")

ps.its.inall@phy_tree <- tree
ps.its.inall  <- prune_taxa(taxa_sums(ps.its.inall) > 0, ps.its.inall)


ps.its.exl  <- phyloseq::filter_taxa(ps_its, function(x) sum(x > 10) < (0.1*length(x)), TRUE)
ps.its.exl  <- prune_taxa(taxa_sums(ps.its.exl) > 100, ps.its.exl)
ps.its.exl.taxa <- taxa_names(ps.its.exl)

amp_its_inall <- phyloseq_to_ampvis2(ps.its.inall)

Base statistics for its

amp_its

## ampvis2 object with 4 elements. 
## Summary of OTU table:
##      Samples         OTUs  Total#Reads    Min#Reads    Max#Reads Median#Reads 
##           36         1295       815134        13533        34994        22657 
##    Avg#Reads 
##     22642.61 
## 
## Assigned taxonomy:
##     Kingdom      Phylum       Class       Order      Family       Genus 
##    554(43%) 453(34.98%) 362(27.95%) 335(25.87%) 321(24.79%) 306(23.63%) 
##     Species 
## 207(15.98%) 
## 
## Metadata variables: 4 
##  SampleID, Day, Description, Group

Major ITS phylotypes

amp_heatmap(amp_its,
            tax_show = 40,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=TRUE,
            showRemainingTaxa = TRUE)

## Warning: Transformation introduced infinite values in discrete y-axis

All by the phylums

amp_heatmap(amp_its,
            tax_show = 7,
            group_by = "Day",
            tax_aggregate = "Phylum",
            tax_add = "Kingdom",
            normalise=TRUE,
            showRemainingTaxa = TRUE)

## Warning: Transformation introduced infinite values in discrete y-axis

There are going to WGCNA pipeline(116 phylotypes)

amp_heatmap(amp_its_inall,
            tax_show = 40 ,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=TRUE,
            showRemainingTaxa = TRUE)

## Warning: Transformation introduced infinite values in discrete y-axis

beta_custom_norm_NMDS_elli_w(ps_its, Group="Day", Color="Day")

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1834317 
## Run 1 stress 0.1834317 
## ... Procrustes: rmse 0.0000201183  max resid 0.00005501659 
## ... Similar to previous best
## Run 2 stress 0.1841584 
## Run 3 stress 0.1834317 
## ... New best solution
## ... Procrustes: rmse 0.000009091179  max resid 0.00003505055 
## ... Similar to previous best
## Run 4 stress 0.1841584 
## Run 5 stress 0.1834317 
## ... New best solution
## ... Procrustes: rmse 0.00002509501  max resid 0.00007022262 
## ... Similar to previous best
## Run 6 stress 0.1834317 
## ... Procrustes: rmse 0.000003875879  max resid 0.00001214466 
## ... Similar to previous best
## Run 7 stress 0.1834317 
## ... Procrustes: rmse 0.000006013033  max resid 0.00002521067 
## ... Similar to previous best
## Run 8 stress 0.1841585 
## Run 9 stress 0.1834317 
## ... Procrustes: rmse 0.000005395813  max resid 0.00001721888 
## ... Similar to previous best
## Run 10 stress 0.1841584 
## Run 11 stress 0.2180006 
## Run 12 stress 0.1834317 
## ... Procrustes: rmse 0.0000189146  max resid 0.00008526694 
## ... Similar to previous best
## Run 13 stress 0.2171479 
## Run 14 stress 0.1841584 
## Run 15 stress 0.2097497 
## Run 16 stress 0.1834317 
## ... Procrustes: rmse 0.0000264768  max resid 0.00007101285 
## ... Similar to previous best
## Run 17 stress 0.1893642 
## Run 18 stress 0.1834317 
## ... Procrustes: rmse 0.000009293648  max resid 0.00002426741 
## ... Similar to previous best
## Run 19 stress 0.1834317 
## ... Procrustes: rmse 0.000008230916  max resid 0.00002983058 
## ... Similar to previous best
## Run 20 stress 0.1834317 
## ... Procrustes: rmse 0.00002722606  max resid 0.00008598357 
## ... Similar to previous best
## *** Solution reached

The first plot is clr normalization (compositional), the second is vst(DESeq2). In contrast to bacteria - for vst WGCNA did not show good results (no plateau yield).
We still applied vst stabilization, but we have to keep in mind that the resulting output is not great.

otu.inall <- as.data.frame(ps.its.inall@otu_table@.Data) 
ps.inall.clr <- ps.its.inall
otu_table(ps.inall.clr) <- phyloseq::otu_table(compositions::clr(otu.inall), taxa_are_rows = FALSE)
data <- ps.inall.clr@otu_table@.Data %>% 
  as.data.frame()
rownames(data) <- as.character(ps.inall.clr@sam_data$Description)
powers <-  c(c(1:10), seq(from = 12, to=30, by=1))
sft <-  pickSoftThreshold(data, powerVector = powers, verbose = 5, networkType = "signed")

## pickSoftThreshold: will use block size 116.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 116 of 116
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.3860 -25.70         0.5550 57.2000 57.100000 59.600
## 2      2   0.4580 -17.60         0.3190 29.5000 29.300000 32.400
## 3      3   0.4860 -11.60         0.3420 15.8000 15.600000 18.500
## 4      4   0.6370  -7.48         0.5340  8.7100  8.530000 11.300
## 5      5   0.7900  -5.30         0.7300  4.9900  4.780000  7.370
## 6      6   0.8050  -3.89         0.7500  2.9600  2.770000  5.090
## 7      7   0.8400  -2.94         0.8030  1.8300  1.650000  3.710
## 8      8   0.7890  -2.65         0.7530  1.1700  1.000000  2.820
## 9      9   0.2080  -7.17        -0.0165  0.7790  0.626000  2.230
## 10    10   0.8000  -1.88         0.8250  0.5400  0.401000  1.810
## 11    12   0.8610  -1.51         0.9500  0.2890  0.173000  1.360
## 12    13   0.7990  -1.44         0.8150  0.2230  0.121000  1.220
## 13    14   0.8000  -1.34         0.8260  0.1770  0.083400  1.120
## 14    15   0.7480  -1.32         0.7690  0.1440  0.057500  1.020
## 15    16   0.1540  -2.22        -0.0813  0.1190  0.040100  0.941
## 16    17   0.1530  -2.10        -0.0816  0.1010  0.027800  0.869
## 17    18   0.8570  -1.15         0.8780  0.0866  0.019400  0.805
## 18    19   0.0897  -1.46        -0.0687  0.0754  0.013700  0.746
## 19    20   0.0930  -1.43        -0.0720  0.0663  0.009910  0.693
## 20    21   0.2140  -2.01         0.0509  0.0589  0.007180  0.660
## 21    22   0.2200  -2.04         0.0648  0.0527  0.005230  0.643
## 22    23   0.0951  -1.76        -0.0800  0.0475  0.003800  0.627
## 23    24   0.1930  -1.80        -0.0248  0.0431  0.002730  0.612
## 24    25   0.1960  -1.82        -0.0191  0.0393  0.001970  0.598
## 25    26   0.2010  -2.33        -0.0126  0.0360  0.001440  0.585
## 26    27   0.2220  -2.37         0.0609  0.0331  0.001070  0.572
## 27    28   0.2290  -2.40         0.0789  0.0305  0.000788  0.559
## 28    29   0.2320  -2.67         0.0387  0.0283  0.000584  0.547
## 29    30   0.2430  -2.73         0.0567  0.0263  0.000433  0.536

plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")

ps.varstab <- ps_vst(ps.its.inall, "Day")

data2 <- ps.varstab@otu_table@.Data %>% 
  as.data.frame()

rownames(data2) <- as.character(ps.varstab@sam_data$Description)
powers <-  c(seq(from = 1, to=10, by=0.5), seq(from = 11, to=20, by=1))
sft2 <-  pickSoftThreshold(data2, powerVector = powers, verbose = 5, networkType = "signed hybrid")

## pickSoftThreshold: will use block size 116.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 116 of 116
##    Power SFT.R.sq  slope truncated.R.sq mean.k.  median.k. max.k.
## 1    1.0    0.352 -1.180        0.82400 11.7000 11.1000000  21.40
## 2    1.5    0.511 -0.965        0.80800  6.4500  5.9700000  13.20
## 3    2.0    0.688 -1.020        0.82800  3.9000  3.3000000   9.04
## 4    2.5    0.787 -0.915        0.80900  2.5200  1.9300000   6.41
## 5    3.0    0.828 -0.797        0.78200  1.7200  1.1700000   4.69
## 6    3.5    0.807 -0.913        0.75900  1.2300  0.7570000   3.81
## 7    4.0    0.789 -1.110        0.73200  0.9160  0.5010000   3.54
## 8    4.5    0.200 -2.790        0.08170  0.7060  0.3370000   3.34
## 9    5.0    0.218 -2.720        0.11300  0.5610  0.2300000   3.18
## 10   5.5    0.235 -3.520        0.08280  0.4570  0.1590000   3.04
## 11   6.0    0.165 -2.140       -0.07370  0.3800  0.1100000   2.92
## 12   6.5    0.204 -2.250       -0.00278  0.3220  0.0758000   2.82
## 13   7.0    0.208 -2.210        0.00364  0.2780  0.0538000   2.73
## 14   7.5    0.192 -2.760       -0.03940  0.2430  0.0391000   2.65
## 15   8.0    0.192 -2.670       -0.03910  0.2150  0.0280000   2.57
## 16   8.5    0.196 -2.600       -0.03270  0.1920  0.0203000   2.50
## 17   9.0    0.228 -2.770        0.01890  0.1740  0.0148000   2.44
## 18   9.5    0.213 -2.990       -0.00815  0.1580  0.0108000   2.37
## 19  10.0    0.214 -3.060       -0.00549  0.1450  0.0078000   2.32
## 20  11.0    0.230 -3.170        0.01010  0.1240  0.0041500   2.21
## 21  12.0    0.262 -3.210        0.05580  0.1080  0.0024100   2.12
## 22  13.0    0.222 -3.130        0.01010  0.0955  0.0014200   2.03
## 23  14.0    0.234 -3.060        0.02160  0.0855  0.0008090   1.95
## 24  15.0    0.259 -3.140        0.04770  0.0774  0.0004490   1.88
## 25  16.0    0.270 -3.070        0.06090  0.0707  0.0002500   1.82
## 26  17.0    0.233 -2.840        0.04610  0.0651  0.0001400   1.76
## 27  18.0    0.244 -2.790        0.05240  0.0602  0.0000789   1.71
## 28  19.0    0.251 -2.860        0.06040  0.0561  0.0000444   1.66
## 29  20.0    0.280 -3.030        0.08580  0.0525  0.0000250   1.62

plot(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")

# WGCNA, as old and not supported
detachAllPackages()
library(WGCNA)

net3 <- WGCNA::blockwiseModules(data2,
                          power=3,
                          TOMType="signed",
                          networkType="signed hybrid",
                          nThreads=0)


mergedColors2 <-  WGCNA::labels2colors(net3$colors, colorSeq = c("salmon", "darkgreen", "cyan", "red", "blue", "plum"))

plotDendroAndColors(
  net3$dendrograms[[1]],
  mergedColors2[net3$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE,
  hang = 0.03,
  addGuide = TRUE,
  guideHang = 0.05)

library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(heatmaply)
library(WGCNA)
library(phyloseq)
library(ggtree)
library(tidyverse)
library(KneeArrower)

The green cluster is the background, which in the case of the red cluster for bacteria that here, I don’t think it can be called a cluster at all.
WGCNA as an analysis helps clustering based on correlations by finding a soft threshold(power, see plot where the red numbers still can’t reach the plateau),
the darkgreen and red cluster are simply unconnected phylotypes.

modules_of_interest = mergedColors2 %>% 
  unique()

module_df <- data.frame(
  asv = names(net3$colors),
  colors = mergedColors2
)

# module_df[module_df == "yellow"] <- "salmon"

submod <-  module_df %>%
  subset(colors %in% modules_of_interest)

row.names(module_df) = module_df$asv

subexpr = as.data.frame(t(data2))[submod$asv,]

submod_df <- data.frame(subexpr) %>%
  mutate(
    asv = row.names(.)
  ) %>%
  pivot_longer(-asv) %>%
  mutate(
    module = module_df[asv,]$colors
  )

submod_df <- submod_df %>% 
  mutate(name = gsub("\\_.*","",submod_df$name)) %>% 
  group_by(name, asv) %>% 
  summarise(value = mean(value), asv = asv, module = module) %>% 
  relocate(c(asv, name, value, module)) %>% 
  ungroup() %>% 
  mutate(module = as.factor(module))

p <- submod_df %>% 
  ggplot(., aes(x=name, y=value, group=asv)) +
  geom_line(aes(color = module),
            alpha = 0.2) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.position = "none") +
  facet_grid(rows = vars(module)) +
  labs(x = "day",
       y = "normalized abundance")

p + scale_color_manual(values = levels(submod_df$module)) +
  scale_x_discrete(labels=(c("D01"="3",
                                "D03"="14",
                                "D05"="28",
                                "D07"="49" ,
                                "D08"="63",
                                "D10"="91",
                                "D12"="119",
                                "D13"="140",
                                "D14"="161",
                                "D15"="182")
         ))

l_its <- color_filt_broken(ps_its, submod_df, ps.its.inall)
l_its

$cyan \(cyan\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 36 taxa and 36 samples ] sample_data() Sample Data: [ 36 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 36 taxa by 7 taxonomic ranks ] refseq() DNAStringSet: [ 36 reference sequences ]

\(cyan\)amp ampvis2 object with 4 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 36 36 318660 0 23584 9007.5 Avg#Reads 8851.67

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 28(78%) 25(69.44%) 17(47.22%) 12(33.33%) 12(33.33%) 12(33.33%) 9(25%)

Metadata variables: 4 SampleID, Day, Description, Group

\(cyan\)heat \(cyan\)heat_rel \(cyan\)tree \(cyan\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq1	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Chaetosphaeriales	f__Chaetosphaeriaceae	g__Chloridium	s__aseptatum
Seq2	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	s__verticillata
Seq12	k__Fungi	p__Ascomycota	c__Sordariomycetes	NA	NA	NA	NA
Seq184	NA	NA	NA	NA	NA	NA	NA
Seq10	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	s__canina
Seq73	k__Alveolata	NA	NA	NA	NA	NA	NA
Seq5	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	s__verticillata
Seq118	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq18	k__Fungi	p__Ascomycota	c__Sordariomycetes	NA	NA	NA	NA
Seq93	k__Fungi	p__Ascomycota	c__Dothideomycetes	o__Venturiales	f__Sympoventuriaceae	g__Scolecobasidium	s__constrictum
Seq78	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq101	NA	NA	NA	NA	NA	NA	NA
Seq4	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	s__verticillata
Seq62	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq48	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq36	k__Fungi	p__Ascomycota	c__Sordariomycetes	NA	NA	NA	NA
Seq44	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq103	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq166	k__Fungi	p__Ascomycota	c__Eurotiomycetes	o__Chaetothyriales	f__Herpotrichiellaceae	g__Exophiala	NA
Seq231	k__Eukaryota_kgd_Incertae_sedis	NA	NA	NA	NA	NA	NA
Seq21	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Clavicipitaceae	g__Metarhizium	s__marquandii
Seq108	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq74	k__Fungi	p__Ascomycota	c__Sordariomycetes	NA	NA	NA	NA
Seq195	NA	NA	NA	NA	NA	NA	NA
Seq28	k__Fungi	p__Ascomycota	c__Leotiomycetes	o__Helotiales	f__Helotiaceae	g__Scytalidium	NA
Seq104	k__Metazoa	p__Nematoda	c__Chromadorea	o__Rhabditida	f__Cephalobidae	g__Pseudacrobeles	NA
Seq29	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Glomerellales	f__Glomerellaceae	g__Colletotrichum	s__sidae
Seq46	k__Fungi	p__Ascomycota	c__Sordariomycetes	NA	NA	NA	NA
Seq122	NA	NA	NA	NA	NA	NA	NA
Seq235	NA	NA	NA	NA	NA	NA	NA
Seq191	NA	NA	NA	NA	NA	NA	NA
Seq94	k__Alveolata	NA	NA	NA	NA	NA	NA
Seq76	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Sordariales_fam_Incertae_sedis	g__Staphylotrichum	s__boninense
Seq155	NA	NA	NA	NA	NA	NA	NA
Seq190	NA	NA	NA	NA	NA	NA	NA
Seq32	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA

$darkgreen \(darkgreen\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 56 taxa and 36 samples ] sample_data() Sample Data: [ 36 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 56 taxa by 7 taxonomic ranks ] refseq() DNAStringSet: [ 56 reference sequences ]

\(darkgreen\)amp ampvis2 object with 4 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 36 56 261584 1790 20315 6577 Avg#Reads 7266.22

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 52(93%) 51(91.07%) 48(85.71%) 48(85.71%) 48(85.71%) 48(85.71%) 34(60.71%)

Metadata variables: 4 SampleID, Day, Description, Group

\(darkgreen\)heat \(darkgreen\)heat_rel \(darkgreen\)tree \(darkgreen\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq92	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Bolbitiaceae	g__Conocybe	s__zeylanica
Seq49	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Bolbitiaceae	g__Conocybe	s__zeylanica
Seq91	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	s__canina
Seq6	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Chaetomiaceae	g__Humicola	s__sardiniae
Seq16	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Bolbitiaceae	g__Conocybe	s__zeylanica
Seq67	k__Metazoa	p__Annelida	c__Clitellata	o__Enchytraeida	f__Enchytraeidae	g__Fridericia	NA
Seq13	k__Fungi	p__Ascomycota	c__Eurotiomycetes	o__Eurotiales	f__Trichocomaceae	g__Talaromyces	NA
Seq7	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Bolbitiaceae	g__Conocybe	s__zeylanica
Seq25	k__Fungi	p__Basidiomycota	c__Cystobasidiomycetes	o__Cystobasidiales	f__Cystobasidiaceae	g__Occultifur	NA
Seq69	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Nectriaceae	g__Gibberella	s__intricans
Seq86	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Bolbitiaceae	g__Conocybe	s__zeylanica
Seq205	k__Fungi	p__Mucoromycota	c__Mucoromycetes	o__Mucorales	f__Rhizopodaceae	g__Rhizopus	s__arrhizus
Seq50	k__Fungi	p__Basidiomycota	c__Cystobasidiomycetes	o__Cystobasidiales	f__Cystobasidiaceae	g__Occultifur	NA
Seq193	k__Heterolobosa	p__Heterolobosa_phy_Incertae_sedis	c__Heterolobosea	o__Schizopyrenida	f__Vahlkampfiidae	g__Naegleria	NA
Seq136	k__Fungi	p__Ascomycota	c__Eurotiomycetes	o__Chaetothyriales	f__Herpotrichiellaceae	g__Exophiala	NA
Seq23	k__Fungi	p__Mucoromycota	c__Mucoromycetes	o__Mucorales	f__Mucoraceae	g__Actinomucor	s__elegans
Seq175	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Bionectriaceae	g__Clonostachys	s__rosea
Seq281	k__Heterolobosa	p__Heterolobosa_phy_Incertae_sedis	c__Heterolobosea	o__Schizopyrenida	f__Vahlkampfiidae	g__Allovahlkampfia	NA
Seq47	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Nectriaceae	g__Gibberella	s__intricans
Seq137	k__Fungi	p__Ascomycota	c__Dothideomycetes	o__Venturiales	f__Sympoventuriaceae	g__Ochroconis	s__tshawytschae
Seq9	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Stachybotryaceae	g__Albifimbria	s__verrucaria
Seq42	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Bolbitiaceae	g__Conocybe	s__zeylanica
Seq45	k__Fungi	NA	NA	NA	NA	NA	NA
Seq11	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Nectriaceae	g__Gibberella	s__intricans
Seq89	k__Fungi	p__Basidiomycota	c__Cystobasidiomycetes	o__Cystobasidiales	f__Cystobasidiaceae	g__Occultifur	NA
Seq171	k__Fungi	p__Ascomycota	c__Dothideomycetes	o__Venturiales	f__Sympoventuriaceae	g__Ochroconis	s__tshawytschae
Seq178	k__Fungi	p__Ascomycota	c__Dothideomycetes	o__Venturiales	f__Sympoventuriaceae	g__Scolecobasidium	s__constrictum
Seq22	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	s__verticillata
Seq31	k__Fungi	p__Mucoromycota	c__Mucoromycetes	o__Mucorales	f__Mucoraceae	g__Actinomucor	s__elegans
Seq65	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	NA
Seq75	NA	NA	NA	NA	NA	NA	NA
Seq85	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq54	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	NA
Seq121	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Nectriaceae	g__Gibberella	s__intricans
Seq125	k__Fungi	p__Ascomycota	c__Dothideomycetes	o__Venturiales	f__Sympoventuriaceae	g__Ochroconis	s__tshawytschae
Seq181	k__Fungi	p__Ascomycota	c__Eurotiomycetes	o__Eurotiales	f__Aspergillaceae	g__Penicillium	s__bialowiezense
Seq24	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Nectriaceae	g__Gibberella	s__intricans
Seq27	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Stachybotryaceae	g__Stachybotrys	s__chartarum
Seq66	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Stachybotryaceae	g__Albifimbria	s__verrucaria
Seq114	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Bolbitiaceae	g__Conocybe	s__zeylanica
Seq15	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Nectriaceae	g__Gibberella	s__intricans
Seq226	k__Fungi	p__Ascomycota	c__Dothideomycetes	o__Venturiales	f__Sympoventuriaceae	g__Scolecobasidium	s__constrictum
Seq26	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Chaetomiaceae	g__Chaetomium	s__iranianum
Seq59	k__Fungi	p__Ascomycota	c__Dothideomycetes	o__Venturiales	f__Sympoventuriaceae	g__Ochroconis	s__tshawytschae
Seq41	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Cantharellales	f__Ceratobasidiaceae	g__Waitea	s__circinata
Seq157	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Chaetomiaceae	g__Zopfiella	NA
Seq124	k__Fungi	p__Ascomycota	c__Eurotiomycetes	o__Eurotiales	f__Aspergillaceae	g__Aspergillus	NA
Seq254	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq17	k__Metazoa	p__Nematoda	NA	NA	NA	NA	NA
Seq39	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Coniochaetales	f__Coniochaetaceae	g__Coniochaeta	s__verticillata
Seq180	NA	NA	NA	NA	NA	NA	NA
Seq20	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Chaetomiaceae	g__Chaetomium	NA
Seq102	NA	NA	NA	NA	NA	NA	NA
Seq70	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Sordariales_fam_Incertae_sedis	g__Staphylotrichum	s__boninense
Seq150	k__Heterolobosa	p__Heterolobosa_phy_Incertae_sedis	c__Heterolobosea	o__Schizopyrenida	f__Vahlkampfiidae	g__Allovahlkampfia	NA
Seq132	NA	NA	NA	NA	NA	NA	NA

$salmon \(salmon\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 24 taxa and 36 samples ] sample_data() Sample Data: [ 36 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 24 taxa by 7 taxonomic ranks ] refseq() DNAStringSet: [ 24 reference sequences ]

\(salmon\)amp ampvis2 object with 4 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 36 24 121335 20 14955 1772 Avg#Reads 3370.42

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 22(92%) 18(75%) 15(62.5%) 14(58.33%) 13(54.17%) 13(54.17%) 11(45.83%)

Metadata variables: 4 SampleID, Day, Description, Group

\(salmon\)heat \(salmon\)heat_rel \(salmon\)tree \(salmon\)taxa

	Kingdom	Phylum	Class	Order	Family	Genus	Species
Seq3	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Lasiosphaeriaceae	g__Schizothecium	s__inaequale
Seq19	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Psathyrellaceae	g__Coprinellus	s__flocculosus
Seq53	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Stachybotryaceae	g__Albifimbria	s__verrucaria
Seq99	k__Viridiplantae	p__Anthophyta	NA	NA	NA	NA	NA
Seq215	k__Viridiplantae	NA	NA	NA	NA	NA	NA
Seq30	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Chaetomiaceae	g__Humicola	s__sardiniae
Seq43	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Lasiosphaeriaceae	g__Schizothecium	NA
Seq68	NA	NA	NA	NA	NA	NA	NA
Seq63	k__Fungi	p__Mucoromycota	c__Mucoromycetes	o__Mucorales	f__Mucoraceae	g__Actinomucor	s__elegans
Seq100	k__Fungi	NA	NA	NA	NA	NA	NA
Seq140	k__Fungi	p__Ascomycota	c__Dothideomycetes	o__Pleosporales	NA	NA	NA
Seq33	k__Fungi	p__Basidiomycota	c__Agaricomycetes	o__Agaricales	f__Psathyrellaceae	g__Coprinellus	s__flocculosus
Seq245	NA	NA	NA	NA	NA	NA	NA
Seq138	k__Viridiplantae	p__Anthophyta	c__Eudicotyledonae	NA	NA	NA	NA
Seq8	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Lasiosphaeriaceae	g__Schizothecium	s__inaequale
Seq79	k__Viridiplantae	p__Anthophyta	NA	NA	NA	NA	NA
Seq161	k__Eukaryota_kgd_Incertae_sedis	NA	NA	NA	NA	NA	NA
Seq96	k__Metazoa	p__Nematoda	c__Chromadorea	o__Rhabditida	f__Cephalobidae	g__Acrobeloides	s__nanus
Seq147	k__Fungi	NA	NA	NA	NA	NA	NA
Seq123	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Chaetomiaceae	g__Chaetomium	s__jodhpurense
Seq152	k__Viridiplantae	p__Anthophyta	NA	NA	NA	NA	NA
Seq52	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Sordariales	f__Lasiosphaeriaceae	g__Schizothecium	NA
Seq40	k__Fungi	p__Mucoromycota	c__Mucoromycetes	o__Mucorales	f__Mucoraceae	g__Actinomucor	s__elegans
Seq37	k__Fungi	p__Ascomycota	c__Sordariomycetes	o__Hypocreales	f__Nectriaceae	g__Gibberella	s__intricans

p.observed <- plot_alpha_w_toc(ps = ps_its, group = "Group", metric = c("Observed")) + 
  theme(axis.title.y = element_blank())

p.shannon <- plot_alpha_w_toc(ps = ps_its, group = "Group", metric = c("Shannon")) + 
  theme(axis.title.y = element_blank())

p.simpson <- plot_alpha_w_toc(ps = ps_its, group = "Group", metric = c("InvSimpson")) + 
  theme(axis.title.y = element_blank())

ggpubr::ggarrange(p.observed, p.shannon, p.simpson, ncol = 3)

p.observed <- plot_alpha_w_toc(ps = ps_its, group = "Day", metric = c("Observed")) + 
  theme(axis.title.y = element_blank())

p.shannon <- plot_alpha_w_toc(ps = ps_its, group = "Day", metric = c("Shannon")) + 
  theme(axis.title.y = element_blank())

p.simpson <- plot_alpha_w_toc(ps = ps_its, group = "Day", metric = c("InvSimpson")) + 
  theme(axis.title.y = element_blank())

ggpubr::ggarrange(p.observed, p.shannon, p.simpson, ncol = 3)

ps.m <- phyloseq::psmelt(ps_its)
ps.m <- ps.m %>% 
  mutate_if(is.character, as.factor) 

ps.data.out <- ps.m %>%
  select(-Group) %>% 
  pivot_wider(names_from = c(Day, Description, Sample), values_from = Abundance, values_fill = 0) 

#create empty dataframe with columnnames
external_empty_dataframe <- data.frame(OTU=factor(), cluster=factor(), stringsAsFactors = FALSE)

for (i in names(l_its)) {
  a <- taxa_names(l_its[[i]][["ps"]])
  b <- rep(i, length(a))
  d <- data.frame(OTU = as.factor(a),
                  cluster = as.factor(b))
  external_empty_dataframe <- rbind(external_empty_dataframe, d)
}

clusters.otu.df <- external_empty_dataframe
# add exl taxa -- taxa "exl"
d <- data.frame(OTU = as.factor(ps.its.exl.taxa),
                cluster = as.factor(rep("exl", length(ps.its.exl.taxa))))

clusters.otu.df <- rbind(clusters.otu.df, d)
ps.data.out.exl <- left_join(clusters.otu.df, ps.data.out, by="OTU")
# write.table(ps.data.out.exl, file = "ps.its.data.out.tsv", sep = "\t")

0.11 Packages info

sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3         reshape2_1.4.4        KneeArrower_1.0.0    
##  [4] ggtree_3.4.2          heatmaply_1.3.0       viridis_0.6.2        
##  [7] viridisLite_0.4.1     plotly_4.10.0         ampvis2_2.7.28       
## [10] ggpubr_0.4.0          forcats_0.5.2         stringr_1.4.1        
## [13] dplyr_1.0.9           purrr_0.3.4           readr_2.1.2          
## [16] tidyr_1.2.0           tibble_3.1.8          ggplot2_3.3.6        
## [19] tidyverse_1.3.2       phyloseq_1.40.0       WGCNA_1.71           
## [22] fastcluster_1.2.3     dynamicTreeCut_1.63-1
## 
## loaded via a namespace (and not attached):
##   [1] bit64_4.0.5                 knitr_1.39                 
##   [3] DelayedArray_0.22.0         data.table_1.14.2          
##   [5] rpart_4.1.16                KEGGREST_1.36.3            
##   [7] RCurl_1.98-1.8              doParallel_1.0.17          
##   [9] generics_0.1.3              BiocGenerics_0.42.0        
##  [11] preprocessCore_1.58.0       cowplot_1.1.1              
##  [13] microbiome_1.18.0           RSQLite_2.2.16             
##  [15] proxy_0.4-27                bit_4.0.4                  
##  [17] tzdb_0.3.0                  bayesm_3.1-4               
##  [19] webshot_0.5.3               xml2_1.3.3                 
##  [21] lubridate_1.8.0             SummarizedExperiment_1.26.1
##  [23] assertthat_0.2.1            gargle_1.2.0               
##  [25] xfun_0.32                   hms_1.1.2                  
##  [27] jquerylib_0.1.4             evaluate_0.16              
##  [29] TSP_1.2-1                   DEoptimR_1.0-11            
##  [31] fansi_1.0.3                 dendextend_1.16.0          
##  [33] dbplyr_2.2.1                readxl_1.4.1               
##  [35] igraph_1.3.4                DBI_1.1.3                  
##  [37] geneplotter_1.74.0          htmlwidgets_1.5.4          
##  [39] tensorA_0.36.2              googledrive_2.0.0          
##  [41] stats4_4.2.0                ellipsis_0.3.2             
##  [43] corrplot_0.92               backports_1.4.1            
##  [45] energy_1.7-10               signal_0.7-7               
##  [47] permute_0.9-7               picante_1.8.2              
##  [49] annotate_1.74.0             compositions_2.0-4         
##  [51] deldir_1.0-6                MatrixGenerics_1.8.1       
##  [53] vctrs_0.4.1                 Biobase_2.56.0             
##  [55] abind_1.4-5                 cachem_1.0.6               
##  [57] withr_2.5.0                 ggforce_0.3.4              
##  [59] robustbase_0.95-0           checkmate_2.1.0            
##  [61] treeio_1.20.2               vegan_2.6-2                
##  [63] mnormt_2.1.0                cluster_2.1.3              
##  [65] gsl_2.1-7.1                 ape_5.6-2                  
##  [67] lazyeval_0.2.2              crayon_1.5.1               
##  [69] genefilter_1.78.0           pkgconfig_2.0.3            
##  [71] labeling_0.4.2              tweenr_2.0.1               
##  [73] GenomeInfoDb_1.32.3         nlme_3.1-157               
##  [75] seriation_1.3.6             nnet_7.3-17                
##  [77] rlang_1.0.4                 lifecycle_1.0.1            
##  [79] registry_0.5-1              modelr_0.1.9               
##  [81] cellranger_1.1.0            polyclip_1.10-0            
##  [83] matrixStats_0.62.0          rngtools_1.5.2             
##  [85] Matrix_1.4-1                aplot_0.1.6                
##  [87] carData_3.0-5               Rhdf5lib_1.18.2            
##  [89] boot_1.3-28                 reprex_2.0.2               
##  [91] base64enc_0.1-3             googlesheets4_1.0.1        
##  [93] png_0.1-7                   rootSolve_1.8.2.3          
##  [95] bitops_1.0-7                rhdf5filters_1.8.0         
##  [97] Biostrings_2.64.1           blob_1.2.3                 
##  [99] doRNG_1.8.2                 gridGraphics_0.5-1         
## [101] jpeg_0.1-9                  rstatix_0.7.0              
## [103] S4Vectors_0.34.0            ggsignif_0.6.3             
## [105] scales_1.2.1                memoise_2.0.1              
## [107] magrittr_2.0.3              plyr_1.8.7                 
## [109] zlibbioc_1.42.0             compiler_4.2.0             
## [111] RColorBrewer_1.1-3          DESeq2_1.36.0              
## [113] cli_3.3.0                   ade4_1.7-19                
## [115] ANCOMBC_1.6.2               XVector_0.36.0             
## [117] patchwork_1.1.2             htmlTable_2.4.1            
## [119] Formula_1.2-4               MASS_7.3-57                
## [121] mgcv_1.8-40                 tidyselect_1.1.2           
## [123] stringi_1.7.8               highr_0.9                  
## [125] yaml_2.3.5                  locfit_1.5-9.6             
## [127] latticeExtra_0.6-30         ggrepel_0.9.1              
## [129] grid_4.2.0                  sass_0.4.2                 
## [131] tools_4.2.0                 lmom_2.9                   
## [133] parallel_4.2.0              rstudioapi_0.14            
## [135] foreach_1.5.2               foreign_0.8-82             
## [137] gld_2.6.5                   farver_2.1.1               
## [139] Rtsne_0.16                  digest_0.6.29              
## [141] Rcpp_1.0.9                  GenomicRanges_1.48.0       
## [143] car_3.1-0                   broom_1.0.0                
## [145] httr_1.4.4                  AnnotationDbi_1.58.0       
## [147] psych_2.2.5                 Rdpack_2.4                 
## [149] colorspace_2.0-3            rvest_1.0.3                
## [151] XML_3.99-0.10               fs_1.5.2                   
## [153] IRanges_2.30.1              splines_4.2.0              
## [155] yulab.utils_0.0.5           tidytree_0.4.0             
## [157] expm_0.999-6                multtest_2.52.0            
## [159] Exact_3.1                   ggplotify_0.1.0            
## [161] xtable_1.8-4                jsonlite_1.8.0             
## [163] nloptr_2.0.3                ggfun_0.0.6                
## [165] R6_2.5.1                    Hmisc_4.7-1                
## [167] pillar_1.8.1                htmltools_0.5.3            
## [169] glue_1.6.2                  fastmap_1.1.0              
## [171] BiocParallel_1.30.3         class_7.3-20               
## [173] codetools_0.2-18            mvtnorm_1.1-3              
## [175] utf8_1.2.2                  lattice_0.20-45            
## [177] bslib_0.4.0                 DescTools_0.99.45          
## [179] GO.db_3.15.0                interp_1.1-3               
## [181] survival_3.3-1              rmarkdown_2.15             
## [183] biomformat_1.24.0           munsell_0.5.0              
## [185] e1071_1.7-11                rhdf5_2.40.0               
## [187] GenomeInfoDbData_1.2.8      iterators_1.0.14           
## [189] impute_1.70.0               haven_2.5.1                
## [191] gtable_0.3.0                rbibutils_2.2.9